A hybrid column generation and simulated annealing algorithm for direct aperture optimization

The purpose of this work was to develop a hybrid column generation (CG) and simulated annealing (SA) algorithm for direct aperture optimization (H-DAO) and to show its effectiveness in generating high quality treatment plans for intensity modulated radiation therapy (IMRT) and mixed photon-electron beam radiotherapy (MBRT). The H-DAO overcomes limitations of the CG-DAO with two features improving aperture selection (branch-feature) and enabling aperture shape changes during optimization (SA-feature). The H-DAO algorithm iteratively adds apertures to the plan. At each iteration, a branch is created for each field provided. First, each branch determines the most promising aperture of its assigned field and adds it to a copy of the current apertures. Afterwards, the apertures of each branch undergo an MU-weight optimization followed by an SA-based simultaneous shape and MU-weight optimization and a second MU-weight optimization. The next H-DAO iteration continues the branch with the lowest objective function value. IMRT and MBRT treatment plans for an academic, a brain and a head and neck case generated using the CG-DAO and H-DAO were compared. For every investigated case and both IMRT and MBRT, the H-DAO leads to a faster convergence of the objective function value with number of apertures compared to the CG-DAO. In particular, the H-DAO needs about half the apertures to reach the same objective function value as the CG-DAO. The average aperture areas are 27% smaller for H-DAO than for CG-DAO leading to a slightly larger discrepancy between optimized and final dose. However, a dosimetric benefit remains. The H-DAO was successfully developed and applied to IMRT and MBRT. The faster convergence with number of apertures of the H-DAO compared to the CG-DAO allows to select a better compromise between plan quality and number of apertures.


Introduction
Treatment plans are generated and delivered in photon intensity modulated radiation therapy (IMRT) (Bortfeld 2006) to achieve a highly conformal dose distribution to the target volume. This was enabled through the introduction of the multileaf collimator (MLC) (Convery and Webb 1992), which collimates photon beams delivered from static beam directions in a step-and-shoot or dynamic movement manner such as slidingwindow. Optimization algorithms were introduced to determine a suitable intensity modulation of the photon fields, which are discretized into beamlets (Webb 1989, Bortfeld 2006. These algorithms are classified in two categories: fluence map optimization (FMO) and direct aperture optimization (DAO). FMO optimizes the nonnegative and independent weights of the beamlets in terms of monitor units (MUs) of each provided treatment field simultaneously. FMO is a large-scale convex optimization problem that can be efficiently solved using deterministic algorithms such as gradient descent. However, FMO results only in a non-deliverable fluence map for each field. Thus, a post-processing leaf-sequencing step is added translating the fluence maps to deliverable plans, resulting in a degraded plan quality . In contrast, DAO directly considers the machine constraints of the MLC such as leaf movement constraints or minimal gaps. Thus, DAO deals with MU weighted mechanically deliverable apertures that describe which beamlets are covered by the MLC leaves. However, this leads to a difficult large-scale non-convex optimization problem which cannot be solved efficiently. Several algorithms were developed for DAO: • Column generation based DAO (CG-DAO) (Romeijn et al 2005, Preciado-Walters et al 2006, Men et al 2007, Carlsson 2008, Renaud et al 2017 iteratively adds the most promising aperture shape to the aperture pool of the plan using a pricing mechanism. The most promising aperture shape of all the provided fields is the one with the steepest gradient on the objective function, called price. After each aperture addition, the aperture MU weights are optimized (restricted master problem) using a deterministic algorithm. This approach is computationally efficient and has full freedom in the number of apertures per field. On the other side, the aperture shapes stay fixed as soon as appended to the aperture pool. Moreover, the selection of the most promising aperture is only based on the objective function gradient, i.e. it is unknown by how much the objective function value can be decreased by increasing its MU weight and re-adjusting the MU weights of the already added apertures. Thus, a suboptimal aperture selection is given, especially if many fields are provided as this is the case for 4Pi (Dong et al 2013) or mixed beam radiotherapy (MBRT) (Palma et al 2012, Mueller et al 2017, Renaud et al 2017. • Stochastic DAO approaches such as simulated annealing (SA-DAO) (Shepard et al 2002), quantum tunnel annealing (Pakela et al 2020) and genetic algorithms (Li et al 2003) randomly change the shapes and MU weights of apertures according to a specific scheme. In contrast to the CG-DAO algorithm, the number of apertures per field is pre-defined. Moreover, these approaches start with arbitrary aperture shapes such as conformal to the target or closed. Hence, numerous optimization steps are needed leading to a long computation time. On the other side, these algorithms have basically the ability to overcome local minima, because they also accept changes on the apertures leading to a worse objective function value with a certain probability. In comparison to the CG-DAO, another benefit is the ability to simultaneously optimize the shapes and MU weights of the apertures.
• Local gradient-based leaf refinement approaches such as the direct machine parameter optimization (DMPO) (Hårdemark et al 2004) and aperture shape optimization (ASO) (Cassioli and Unkelbach 2013) typically start with an initial aperture set generated through CG-DAO (Carlsson 2008, Cassioli andUnkelbach 2013) or FMO and leaf-sequencing. Subsequently, they refine the leaf positions locally within the current beamlet using a linear function of the leaf positions approximating the dose distribution. They are able to simultaneously optimize shapes and MU weights of the apertures. However, they can end up in a local minimum making them dependent on the starting conditions.
• Segmentation of fluence maps can also be directly integrated into FMO as shown by Nguyen et al (2017) using a multiphase piecewise constant Mumford-Shah formulation. Thus, the MLC constraints are directly included in the FMO formulation making it a DAO. Potential limitations of the published approach are that it is designed to generate only non-overlapping apertures and that the maximal number of allowed apertures per field is pre-defined.
In this work, a novel algorithm called hybrid DAO (H-DAO) is developed to solve the DAO problem. The H-DAO follows the basic idea of the CG-DAO of adding apertures iteratively, but the H-DAO overcomes limitations of the CG-DAO with the following newly implemented features: • The branch-feature exploring the most promising aperture of each field in a separate branch to identify the aperture improving the objective function value the most.
• The SA-feature applying the simulated annealing algorithm to enable continuous optimization of the aperture shapes and to enable continuous leaf positions not limited by the discrete beamlet grid resolution.
The H-DAO is applied to IMRT and MBRT. We demonstrate the effectiveness of the H-DAO in generating treatment plans of high quality compared to the baseline of a CG-DAO.

Treatment planning process
An implementation of the H-DAO is embedded in the treatment planning process (TPP) as illustrated in figure 1 (top). The TPP considered here is used to create treatment plans for IMRT (photon apertures only) and MBRT (photon and electron apertures) deliverable in a step-and-shoot manner on a TrueBeam (Varian Medical Systems, Palo Alto, CA) treatment unit equipped with a Millennium MLC 120 (Varian Medical Systems, Palo Alto, CA) to collimate photon and electron beams. However, the whole TPP could also be conceptually applied to other MLC based and photon-electron beam supporting treatment units. In the first subprocess, the CT image set is imported into a research version of the Eclipse treatment planning system (TPS) 15.6 (Varian Medical Systems, Palo Alto, CA) and the contouring of the PTV, organs at risk (OARs) and normal tissue (body without PTV) is done using Eclipse.
The second subprocess consists of the manual setup of photon and optionally electron radiation fields within Eclipse. The definition of a field consists of the gantry, table and collimator angles, isocenter location, particle type, beam energy, secondary collimator jaw positions and the beamlet grid size. For photon beams, the secondary collimator jaw positions are equivalent to the maximum rectangular beamlet grid size, while for the While the processes with magenta color are part of the CG-DAO algorithm, the processes with blue and green color are extensions described in this work. (c): Again, the workflow of the H-DAO though presented for a specific optimization example to better illustrate the parallel execution. Each arrow represents a branch run in a separate computer thread. The numbers stand for the number of apertures currently present in the aperture pool. electron beams, the secondary collimator jaws are always set to 15×35 cm 2 due to the electron beam model utilized (Henzen et al 2014). These field definitions stay fixed, but one can setup as many fields as wanted with different field aspects (e.g. for each electron beam direction, all possible beam energies can be setup to give the H-DAO more freedom in selecting the energy of apertures to be added).
To perform Monte Carlo (MC) beamlet dose calculations of the photon and electron fields, the Eclipse interfaced Swiss Monte Carlo Plan (SMCP) (Fix et al 2007) is used. The source of the beamlet dose calculations are pre-simulated phase-spaces located at the treatment head exit plane. The dose distributions of photon and electron beamlets are calculated using Voxel Monte Carlo (VMC++) (Kawrakow and Fippel 2000) and Macro Monte Carlo (MMC) (Neuenschwander and Born 1992, Neuenschwander et al 1995, Fix et al 2013, respectively. Both MC algorithms are embedded within the SMCP framework. The beamlet size is 0.5×0.5 cm 2 or 0.5×1.0 cm 2 depending whether the beamlet belongs to an inner 0.5 cm wide leaf pair or outer thicker 1 cm wide leaf pair. After beamlet dose calculation, the plan optimization is performed using the H-DAO as described in detail in the next section 2.2. Next, a final dose calculation of the optimized apertures considering the impact of the MLC is performed within the SMCP framework. Source of the dose calculation for photon beams is a pre-simulated phase-space located on a plane above the secondary collimator jaws. For photon beams, VMC++ is used to simulate the patient-specific part of the treatment head including secondary collimator jaws and MLC as well as the dose calculation in the patient. For electron beams, the source is a multiple source beam model, called ebm70 After final dose calculation, the MU weights of the apertures are re-optimized using a limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm to reduce the degradation from optimized to final dose distribution, called optimization convergence error (OCE) (Jeraj 2002).

Plan optimization
The goal of the optimization is to minimize the objective function F, which quadratically penalizes deviations of the plan's dose distribution to N DV upper and lower dose-volume objectives ] are the fraction of voxel i overlapping with the considered structure, V , str r , V str s , and V str q , are the summed voxel fractions of the voxels overlapping with the considered structure, M , r M s and M q are the number of voxels of the considered structure, q is the Heaviside function, a r is equal to 1 and −1 for upper and lower dose-volume objectives, respectively, D r is the objected dose and D(V r ) is the dose received by at least the tolerated volume V r of the considered structure for dosevolume objective r, gEUD(t,s) is the gEUD value with tissue-specific factor t, gEUD s is the objected gEUD value, D i q , is the objected normal tissue dose to voxel i which is a function of the nearest distance x i of voxel i to the target, d 0 and ¥ d are the start and end dose parameters, respectively, x start is the start distance of the fall-off and b is the fall-off factor.
Plan dose D i can be formulated as where A ki is the dose delivered to voxel i with unit MU and w k is the MU weight of aperture k. K is the set of all deliverable apertures of the provided fields. For clarification, deliverable does not make any statement about the accuracy of the delivered dose. In this work, deliverable just means mechanically deliverable by the MLC. A ki can be further defined as where F k is the set of all beamlets belonging to the field of aperture k and B ji is the delivered dose to voxel i per unit MU of beamlet j. B ji is calculated by the subprocess 3 of the TPP as described in the previous section 2.1.
] is the transmission factor for beamlet j given by ,, C j L , and C j R , are the position of beamlet j counted from the left and right beamlet grid border, respectively, and P k j L , , and P k j R , , are the left and right leaf position of the leaf pair in aperture k in the line of beamlet j. t Open is 100%, t Tip is 12.9% (photon beam) and 0% (electron beam) and t MLC is 1.3% (photon beam) and 0% (electron beam). These transmission factors of photon beams are an empirically determined model to approximate the transmission through the MLC (t MLC ) and increased transmission through the leaf tip t .
Tip ( ) This model aims to reduce the dose prediction error (DPE) and therefore, the OCE (Jeraj 2002, Bergman et al 2006, Men et al 2007, Mueller et al 2017. In other words, the discrepancy between optimized and final dose distribution is reduced, because the transmission of photon beams through the MLC is already considered during optimization. Beside of that, the transmission model defined above also allows for leaf positions between beamlet borders. The calculation of t G G , Note that equation (7) must be understood as a sum over all possible deliverable apertures. Apertures which are part of the aperture pool have a weight > w 0 k and all other apertures have a weight of 0. The H-DAO algorithm starts with an empty aperture pool, iteratively adds apertures to the pool and optimizes the shapes and MU weights of the apertures at each iteration to minimize F given in equation (1). Figure 1 (middle) illustrates the workflow of the H-DAO algorithm. In the following the details of one iteration are described. One iteration creates as many branches as there are fields provided (branch-feature). The subprocess of these branches are performed in parallel in multiple computer threads. For the current work we used an AMD Epyc2 CPU featuring 64 CPU cores to have enough cores such that no thread needs to share a CPU core with another thread. A specific optimization example of the H-DAO algorithm workflow illustrating the concept of parallel branches is shown in figure 1 (bottom). As visible there, a branch consists of four sequential subprocesses. First, the pricing is performed to find the most promising aperture shape of the field considered in this branch. This determination of the most promising aperture per field is the same as for the CG-DAO (Romeijn et al 2005, Men et al 2007. The most promising aperture is a set of beamlets that can be translated to a deliverable aperture and for which the summed gradient components of the beamlets not covered by the MLC on the objective function, called price, is minimal. The optimization problem to find the most promising aperture of a field can be formulated by For derivation and description of strategies to solve this problem, the work of Romeijn et al (2005) and Men et al (2007) is referenced here.
After determination of the most promising aperture, it is added to a copy of the current aperture pool. It follows an MU weight optimization of the apertures using a projected L-BFGS two loop recursion (Bangert 2011, Nocedal andWright 2013), which is a quasi-Newton algorithm approximating the product of the inverse Hessian and the gradient. The length of the L-BFGS history to store objective function values and gradient values is set to 4. The algorithm is combined with a line search to find an appropriate step length satisfying Armijo's Rule. If any aperture MU weight is below 0 during line search, the MU weight is projected to 0. The L-BFGS terminates if the objective function value is not lowered more than 0.01% in three consecutive iteration.
The next subprocess is the SA-feature, which simultaneously optimizes shapes and weights of the apertures in the pool including the recently added aperture according to an SA cooling schedule. The SA-feature runs maximally for a total of 5000 iterations. It stops earlier, if the objective function value did not decrease more than 0.1% for 250 consecutive iterations. At every iteration, an aperture of the pool is selected randomly for which its shape or MU weight is changed with a probability of P S or (1 − P S ), respectively. The change is accepted if the objective function value is decreased or otherwise with a probability of where T P is the cooling rate, P 0 is the initial value of P, N A is the number of apertures and n S and n W are the number of previous total accepted shape and MU weight changes, respectively. If the shape of the aperture is aimed to be changed, a leaf is randomly selected, and its position is randomly changed according to a normal distribution around the current leaf position and a width of in units of number of beamlets. T S is the cooling rate and σ S0 is the initial width of the normal distribution and N L is the total number of leaf pairs of all apertures in the pool. In case of a weight change, the weight is changed according to a normal distribution around the current aperture weight w k and a width of in relative units of w . k T W is the cooling rate and σ W0 is the initial width of the normal distribution. Following parameter values are used in this work: Note that σ S0 is only 0.3 beamlets wide. Thus, most of the shape changes are small such that the SA-feature follows its purpose of a leaf refinement.
When the SA-feature terminates, the MU weights of the apertures are again optimized using the same L-BFGS implementation. Due to performing an MU weight optimization before and after performing the SAfeature, the SA-feature starts with a better initial solution and the need for more iterations of the SA is also smaller. Note that the MU optimization is convex and involves many optimization variables less (i.e. leaf positions) than the whole DAO optimization problem. Thus, it can be solved very efficiently with the L-BFGS.
When all branches are performed, the objective function value is evaluated for each branch and the aperture pool for the next iteration of the H-DAO is the one with the lowest objective function value among all branches. These H-DAO iterations are repeated until the desired number of apertures in the aperture pool is reached.

Academic and clinical cases
In the computational study of this work, treatment plans of the treatment techniques IMRT and MBRT are generated for an academic case (Mueller et al 2017), a brain case and a head and neck case. Motivation behind the selection of these three cases and the corresponding field setup illustrated in figure 3 is to have an increasing complexity in the geometrical treatment situation and number of fields. The collimator rotation of all the fields is 0°, except the photon fields for the two clinical cases, which have a collimator rotation such that the maximal field width parallel to leaf movement direction is minimized. Reason for this is to get a field width smaller than 15 cm, which is the maximal leaf separation distance between leaves of the same MLC bank.
The academic case is a cylindric homogeneous water phantom with a radius of 10 cm to be treated with a prescribed median dose of 50 Gy to the PTV in 25 fractions. A superficially located PTV including a deep-seated part and two OARs, OAR-lateral and OAR-distal, are contoured. These three structures are extended perpendicular to the transversal plane by 7.4 cm. The first clinical case is a glioblastoma brain case to be treated with a prescribed median dose of 60 Gy to the PTV in 30 fractions. The second clinical case is an oropharynx head and neck case to be treated with a prescribed D95% of 50 Gy to the PTV in 25 fractions. For each case, every optimization used the same dose objectives listed in table 1.

Computational study
A computational study is performed to analyze the influence of the two features added to the CG-DAO algorithm, i.e. the branch-and SA-features, on the following optimization results and properties.  . Photon (top row) and electron (bottom row) field setup for an academic case (left column), a brain case (center column) and a head and neck case (right column) used to create IMRT and MBRT treatment plans. Following structures are visible on the transversal views: PTV (red), OAR-distal (blue) and OAR-lateral (green) for the academic case, PTV (red), brainstem (green), brain (blue), left (cyan) and right (brown) eyes and left (magenta) and right (orange) optic nerves for the brain case, PTV (red), spinal cord (green), right submandibular gland (blue), oral cavity (light brown), larynx (magenta) and pharyngeal constrictors (light yellow) for the head and neck case.
In this computational study, treatment plans of the treatment techniques IMRT and MBRT are generated for the introduced academic and clinical cases using the optimization algorithms listed in table 2. We investigated both IMRT and MBRT to find out if the complexity of the field setup plays any role in the algorithm performance. Due to the choice of two particle types and multiple beam energies, the complexity of the field setup for MBRT can be judged substantially higher. In particular, the fields provided to IMRT are a subset of the fields provided to MBRT, because the same photon fields are provided to both IMRT and MBRT.
To study the convergence with number of apertures of the investigated DAO algorithms, the objective function value after plan optimization (subprocess 4 in the TPP) is collected as a function of the apertures for up to 200 apertures assuming that more than 200 apertures would only lead to marginal improvements in plan Table 1. Objectives used for all the optimizations performed in the computational study. The parameters for the normal-tissue objective have the same values for every case: Start dose d 0 is 95%, end dose ¥ d is 10%, start distance x start of the fall-off is 0.5 cm and fall-off factor b is 0.15. Note that gEUD(t = 1) is equivalent to the mean dose. quality. Furthermore, an FMO is performed with the L-BFGS using the same beamlets without any constraints on smoothness of the fluence map. This FMO optimized plan is used as an ideal benchmark. This is done for IMRT and MBRT for all three cases.
To also see the differences between the different DAO algorithms in the dosimetric space for a specific number of apertures, the TPP is further performed to the end utilizing 50 (academic case), 100 (brain case) and 150 (head and neck case) apertures for all DAO algorithms investigated and the final dose distribution is evaluated by dose-volume-histograms (DVHs) and following dosimetric quantities: PTV dose homogeneity index, mean dose to normal tissue D , mean NT average mean dose to parallel OARs D mean and average D 2% to serial OARs D .

2%
The OARs with a gEUD(t = 1) objective are considered as parallel and the others as serial for the scope of this study. The PTV dose homogeneity index HI is defined as where D p is the prescribed dose and D 2% and D 98% the dose receiving at least 2% and 98% of the PTV volume.
Using the results of these final re-optimized plans, OCE, plan complexity described in numbers of MU and average aperture area A, field contributions (only academic case) and the computation time are also evaluated. The OCE is calculated as where F O is the objective function value after optimization (subprocess 4 of the TPP) and F F is the objective function value after aperture MU weight re-optimization of the final dose distributions (subprocess 6 of the TPP). Connected to the investigation of the OCE, it is evaluated how a potential improvement in terms of objective function value after optimization of the extended DAO algorithms over the CG-DAO evolves after the aperture MU weight re-optimization. Therefore, the improvement after optimization is calculated by A field contribution is defined as the fraction of the PTV mean dose delivered by all the apertures belonging to the corresponding field. In case of MBRT, also photon and electron contributions are analyzed, which is the sum of all photon and electron field contributions, respectively.
To study the statistical uncertainty on the objective function value and the DVH of the DAO algorithms utilizing the SA-feature, the optimizations for the academic case utilizing 50 apertures are repeated 100 times using different seeds to initialize the random number generator. Figure 4 shows the convergence behavior of the objective function value as a function of the number of apertures for all the investigated DAO algorithms. It is visible that for each combination of case and treatment technique (IMRT or MBRT), the fastest convergence is always given by H-DAO, followed by CG-DAO_SA, CG-DAO_Branch and CG-DAO. All the DAO algorithms do not converge to the value given by the FMO, because all the DAO algorithms consider transmission through the MLC in case of photon apertures in contrast to FMO.

Convergence behavior with number of apertures
As the fields provided to IMRT is a subset of the fields provided to MBRT, it is reasonable to hypothesize that the objective function value for MBRT is at least as good as for IMRT. This is true for all the DAO algorithms for at least 10 apertures with the following two exceptions given for the head and neck case: CG-DAO and CG-DAO_SA.
It is noteworthy how many apertures the H-DAO could use less than the CG-DAO to reach the same objective function value, e.g. when the number of apertures of 50, 100 and 150 are selected for the CG-DAO for the academic, brain and head and neck case, respectively, the H-DAO needs on average 44.3% (IMRT) and 55.3% (MBRT) less apertures to reach the same objective function value as the CG-DAO. Figure 5 compares DVHs and table 3 dosimetric quantities of the IMRT and MBRT plans optimized with the different DAO algorithms after performing the whole TPP, i.e. including final dose calculation and MU weight re-optimization. For these plans, specific numbers of apertures are selected, being 50 (academic case), 100 (brain case) and 150 (head and neck case). The DVHs confirm the finding of the convergence behavior analysis that H-DAO performed best followed by CG-DAO_SA, CG-DAO_Branch and CG-DAO. Only the PTV and selected OARs are shown in the DVH comparison for better visibility, but the mentioned finding is in general also confirmed by the DVHs of other structures considered during optimization.

Specific number of apertures
Table 3 also shows that the OCE is higher for the DAO algorithms using the introduced branch-and SAfeature compared to the CG-DAO. However, the OCE is typically small enough such that the improvement after optimization DF o over the CG-DAO is not vanished after aperture weight re-optimization. In case of H-DAO, the improvement DF F is reduced by 15.4% (IMRT) and 14.5% (MBRT) compared to DF O averaged over the three investigated academic and clinical cases.
Regarding plan complexity, table 3 shows that the branch-and SA-features lead to smaller apertures and connected to it, also to higher number of MUs. This is true for every investigated case and both IMRT and MBRT. In case of H-DAO, the aperture areas are 21.8% (IMRT) and 31.2% (MBRT) smaller and the MUs are 11.3% (IMRT) and 30.1% (MBRT) higher compared to CG-DAO averaged over the three investigated academic and clinical cases.
Averaged over the three investigated academic and clinical cases, the CG-DAO_Branch required 19.3% (IMRT) and 27.4% (MBRT) longer computation time than the CG-DAO. Thus, the computation time is not substantially increased due to the use of a multi-core CPU being able to perform the branches in parallel. Again, averaged over the three cases and compared to the CG-DAO, CG-DAO_SA required 146.8% (IMRT) and 110.9% (MBRT) longer computation time and the H-DAO required 179.2% (IMRT) and 201.8% (MBRT) longer computation time. Figure 6 compares the field contributions and photon and electron contribution for the plans generated by the different DAO algorithms utilizing 50 apertures for the academic case. In case of IMRT, the variations between the different DAO algorithms in field contributions are only a few percent. The field contributions also match well with those received by the FMO algorithm. For MBRT, the field contributions vary more between the different DAO algorithms, but they are still within 11% of the FMO field contributions. Overall, the field contributions of the H-DAO match closest with those of the FMO. The same is true for the photon and electron contribution.
The statistical uncertainty of the CG-DAO_SA and H-DAO due to the seed to initialize the random number generator for the SA-feature is demonstrated in figure 7. The standard deviations of the objective function value distributions are 0.0045 (IMRT) and 0.0079 (MBRT) for CG-DAO_SA and 0.0033 (IMRT) and 0.0028 (MBRT) for H-DAO. Thus, the branch-feature included in the H-DAO seems to lead to a lower statistical uncertainty. The dosimetric differences between the worst and best optimization run out of 100 runs are demonstrated in figure 7 for the H-DAO applied to IMRT.

Discussion
The H-DAO algorithm extending the CG-DAO with the branch-and SA-features is successfully implemented and investigated. Each feature alone leads to a faster convergence with number of apertures than the CG-DAO and the combination of the two features also lead to an additional benefit in fast convergence. These statements are generally true for both IMRT and MBRT and for all three investigated cases. The faster convergence with number of apertures can be exploited to create treatment plans with higher dosimetric plan quality or with reduced number of apertures leading to shorter plan delivery time.
Interestingly, CG-DAO and CG-DAO_SA showed a slower convergence with number of apertures for IMRT than for MBRT for the head and neck case up to about 50 apertures, even though the fields provided to IMRT are a subset to those provided to MBRT. This was not observed for CG-DAO_Branch and H-DAO, which both use the branch-feature. Thus, utilizing the branch-feature fulfills better that the provided fields are exploited and not overburdening the DAO algorithm. Therefore, the branch-feature is promising to be also used as a beam angle optimization feature for any CG based DAO algorithm, similarly as it was used for FMO algorithms     Both branch-and SA-feature lead to smaller apertures. Thus, the finding by Cassioli and Unkelbach (2013) that the CG-DAO leads to large apertures, whose MU weight can only be reduced by generating more other apertures, is confirmed. At least for the SA-feature this was expected as it has the purpose of a leaf refinement similar to the aperture shape optimization (ASO) (Cassioli and Unkelbach 2013). Through the evaluation of the OCE, DF o and DF , F it was shown that the smaller apertures can be handled successfully through the whole TPP without substantial deterioration of the objective function value. Main reason for this is a feature considering transmission through the leaves and increased transmission through the leaf ends during optimization. The plans are also deemed to be accurately deliverable by the treatment unit with high accuracy in the delivered dose as this was shown by a recent work of Heath et al about robust optimized MBRT utilizing the H-DAO algorithm for plan generation (Heath et al 2021). There, the delivery of MBRT plans was validated by dose measurements using gafchromic film placed in an anthropomorphic phantom. However, any uncertainties in the modeling of the MLC would be amplified by smaller apertures making an accurate consideration of the treatment head necessary. Furthermore, smaller apertures are usually connected to an increase in MUs and therefore to an increase in delivery time. However, only a minor increase in delivery time is expected when switching from CG-DAO to H-DAO, e.g. for the plans with equal number of apertures shown in this work, the additional delivery time for H-DAO compared to CG-DAO is expected to be about 12 s with 600 MU min −1 due to the additional 120 MUs on average.
The additional computational effort to perform the SA-feature is substantial and even enormous in case of the branch-feature. However, the branch-feature can be easily parallelized by running each branch in a separate thread like this is done in this computational study. Utilizing state-of-the-art CPUs with high number of CPUcores are of high value for this task. The results of this computational study show, that in the worst case the computation time to perform the optimization using H-DAO instead of CG-DAO was not more than 3.9 times increased. Multiple techniques were recently investigated to enhance computational efficiency for DAO that could be combined with the H-DAO.
• Yang et al (2018) replaced the original pricing mechanism with a combination of noise cancellation of the prices using a fuzzy controller followed by aperture generation using threshold segmentation. This allowed to reduce optimization time by 58.61%. This technique could be also used in combination with the H-DAO to determine the most promising aperture per field.
• Men et al (2009) developed a GPU based implementation of the column generation leading to optimization times below 3.8 s. A GPU implementation would be compliant with the branch-feature and the SA-feature could also be implemented to be performed on the GPU.
• MacFarlane et al (2019) reformulated the objective function using a second order Taylor series expansion allowing to find the global minimum by a fast matrix inversion. They applied it to a gradient-based optimization algorithm with at least 70-200 times faster execution and noted that it could be also applied to SA.
• Renaud et al (2017) applied the CG-DAO to MBRT and investigated different aperture addition-schemes adding multiple apertures per iteration. These aperture addition-schemes could be combined with the SAfeature but not directly with the branch-feature presented in this work.
A limitation of this work is the small number of test cases, which does not allow to state any treatment site specific benefits. The H-DAO algorithm could be tuned scenario specific leading potentially to further improvements, e.g. using different parameter values for the SA-feature depending on the treatment site or by adapting the parameter values according to characteristics of the clinical case such as the target size.

Conclusions
The H-DAO algorithm is successfully developed. It extends the CG-DAO algorithm by the branch-feature acting as a more founded decision on the aperture to be added to the aperture pool and by the SA-feature acting as a leaf refinement. This computational study shows that both features lead to a faster convergence of the objective function value with number of apertures. This allows to select a better compromise between dosimetric plan quality and number of apertures.