Comparing variable neighbourhood search algorithms for the direct aperture optimisation in radiotherapy

Introduction

Cancer is a significant global health concern, leading to millions of deaths annually and imposing substantial economic and social burdens. According to GLOBOCAN, there were around $20$ million new cases of cancer in 2022, resulting in nearly $9.7$ million deaths (Bray et al., 2024). This disease ranks as the second leading cause of death worldwide (Ritchie & Roser, 2018), with its prevalence having grown over the past four decades (Quaresma, Coleman & Rachet, 2015). The estimation by the year 2040 considers an increase in diagnosed cancer cases, with an estimated $28.4$ million new cases worldwide, representing a 47% rise from the 2020 figures, assuming a constant rate based on 2020 national estimations (Sung et al., 2021). The escalating incidence of cancer cases is likely to be accompanied by higher mortality rates unless adequate resources are allocated to healthcare systems to manage the growing cancer burden. Primary prevention remains a highly effective approach, as up to half of all cancers are potentially preventable. However, significant efforts are needed to integrate existing effective interventions into healthcare plans (Sung et al., 2021). Recent advancements in research have contributed to improvements in cancer survival rates, which have nearly doubled over the past forty years (Anderson & Nichols, 2020; Quaresma, Coleman & Rachet, 2015).

One of the primary methods for treating cancer is radiotherapy. This treatment involves exposing the patient to high levels of ionising radiation to target and destroy cancer cells. The planning of radiotherapy treatment has evolved significantly, becoming a pillar in the fight against cancer, especially in a global context where this disease represents a considerable burden for society (Bray et al., 2024). In this context, two main approaches stand out: external radiotherapy, which uses machines to focus radiation beams on the tumour from outside the body, and internal radiotherapy or brachytherapy, which involves the direct placement of radioactive material near the tumour for focused dosing, thus minimising the impact on surrounding healthy tissues. Particularly within external radiotherapy, intensity modulated radiation therapy (IMRT) has emerged as a prominent treatment modality, recognised for its ability to precisely direct radiation and tailor it to the shape of the tumour (Daly, 2020). A 2023 report from Emergen Research highlights the growth of the global IMRT market, valued at $2.21 billion in 2022. It is projected to continue expanding, with a compound annual growth rate (CAGR) of 5.4%, driven by factors such as increasing government incentives and preference for non-invasive cancer treatments (Emergen Research, 2023). In this study, we concentrate on IMRT, which administers ionising radiation externally, using a device called a linear accelerator (also called LINAC), from a set of beam angles around the patient’s body in a step-and-shot manner, as depicted in Fig. 1.

IMRT, while clinically effective, encounters a number of challenges. These include the complexity of treatment planning, delivery and quality assurance, which often requires increased time and effort from clinical staff (Ur Rehman et al., 2018). IMRT’s precise dose delivery demands advanced imaging and accurate identification of critical structures, as even slight inaccuracies in defining organ margins can significantly impact treatment due to IMRT’s sensitivity to marginal errors (Beaton et al., 2019). Additionally, IMRT often involves higher costs due to the need for imaging equipment, linear accelerators, synchrotrons, advanced software, computer networks, dosimetry and quality assurance systems, as well as the requirement for a larger, highly trained staff. These challenges emphasise the importance of continuous research and development to enhance the efficiency and precision of IMRT treatment planning.

Given its inherent complexity, the IMRT problem is usually tackled as a set of sequential problems, namely, the beam angle optimisation (BAO), the fluence map optimisation (FMO), and the multi-leaf collimator sequencing (MLC) problems (Cao et al., 2009). The BAO problem looks for the optimal beam angle configuration (BAC). For each identified BAC, the FMO computes the optimal intensities to determine the most effective treatment plan for the corresponding BAC. Ultimately, a collection of achievable aperture configurations for the multi-leaf collimator (MLC) and their associated intensities must be computed (MLC sequencing problem) (Ehrgott et al., 2010).

As mentioned in Salari, Men & Romeijn (2011), one issue with the sequential approach is the large set of apertures generated to solve the MLC sequencing problem. A treatment plan with too many apertures leads to longer beam-on time (BoT), meaning patients are exposed to radiation for extended periods. Longer treatment times are undesirable as they require patients to remain on the couch for extended periods, increasing the risk of inaccuracies in treatment delivery due to patient movement (Dzierma et al., 2014) and reducing the number of patients that can be treated daily. Additionally, the total delivery time of a treatment plan, which comprises the BoT and the time needed by both the linac and the MLC to move between beam angles and achieve the next aperture shape (decomposition time), further contributes to the overall treatment duration (Ahuja & Hamacher, 2005; Baatar et al., 2005). Therefore, it is crucial to decrease the BoT value to shorten the radiation exposure time and reduce the likelihood of treatment interruptions caused by patient discomfort or movement during the procedure.

One strategy treatment planners use to shorten treatment plans is to simplify the treatment plan by rounding the intensities at each beam angle, thereby reducing the number of deliverable aperture shapes required by the treatment plan. However, this rounding strategy also diminishes the quality of the treatment plan (Rocha et al., 2012). As we mentioned before, it is simply not possible to deliver the optimal plan obtained from the FMO problem (as it needs too many aperture shapes). Thus, one strategy proposed in the literature is to include some additional constraints to the FMO problem so that it accounts for the physical considerations of the delivery process. This approach enables us to create a treatment plan that can be directly delivered to the patient without requiring additional steps. The resulting optimisation problem (i.e., the FMO + delivery constraints) is referred to as direct aperture optimisation (DAO) (Shepard et al., 2002).

The DAO problem simultaneously optimises the intensities and shapes of the apertures for each beam in the BAC, ensuring that the treatment plan adheres to the MLC’s physical limitations. By incorporating these constraints during optimisation, DAO eliminates the need for post-optimisation adjustments of the sequential approach, improving the overall quality of the treatment plan (Ludlum & Xia, 2009).

One of the benefits of DAO is its ability to control the number of segments and reduce the BoT required for each treatment plan, thus minimising radiation leakage and the risk of secondary cancers, which is a concern with traditional beamlet-based IMRT methods. According to Broderick, Leech & Coffey (2009), excessive complexity in IMRT plans can increase treatment time and dosimetric uncertainty. Fortunately, DAO can mitigate these issues by simplifying the intensity map and reducing the segments required to achieve high-quality dose distributions.

Aperture shape optimisation and aperture weight optimisation are the two processes involved in DAO (Romeijn et al., 2005). The initial stage involves calculating and adding to the treatment plan the deliverable aperture that offers the greatest potential improvement in the objective function. The aperture weights are usually optimised in the second stage using exact techniques. Following the above approach, DAO produces solutions of significantly higher quality than the sequential approach (Ludlum & Xia, 2009).

The DAO problem has been addressed in various ways. Gradient leaf refining is the foundation of one kind of approach. By determining the link between the objective function and the leaf position and computing its first derivative, these techniques use the leaf position as the optimisation variable. Numerous commercial treatment systems, such as the direct machine parameter optimisation model in the Pinnacle and RayStation systems (Hardemark et al., 2003; Worthy & Wu, 2009), have included these algorithms.

The column generation methods (Romeijn et al., 2005; Preciado-Walters et al., 2006; Carlsson, 2008; Zhang et al., 2019; Salari & Unkelbach, 2013) are another well-known kind of strategies employed in the literature. With this approach, additional deliverable apertures are gradually added to the treatment plan, rather than establishing the initial apertures at the beginning of an iteration. In the iteration phase, the pricing problem is solved to provide a deliverable aperture that can improve the objective function. This aperture is then incorporated into the treatment plan, and the new aperture weights are optimised in the master problem. One aspect to consider is that column generation approaches in the literature converge quickly but do not impose a hard limit on the number of apertures, which could lead to excessively long total treatment times and extremely small apertures (Ripsman et al., 2022).

The DAO problem has also been solved using methods such as stochastic methods (Shepard et al., 2002; Cotrutz & Xing, 2003; Li, Yao & Yao, 2003; Moyano & Cabrera-Guerrero, 2020; Cáceres, Araya & Cabrera-Guerrero, 2021; Fallahi, Mahnam & Niaki, 2021; Fallahi, Mahnam & Akhavan Niaki, 2022; Moyano et al., 2023; Tello-Valenzuela, Moyano & Cabrera-Guerrero, 2023; Tian et al., 2023). These techniques use minor adjustments to the apertures’ leaf positions at each iteration. A modification in leaf position is approved when it enhances the objective function value. One issue with stochastic search and gradient-based leaf refinement techniques is the generation of the initial solution. The quality of the initial solution influences the quality of the given final solution as seen in Pérez Cáceres et al. (2019), Cáceres, Araya & Cabrera-Guerrero (2021), Moyano & Cabrera-Guerrero (2020), Moyano et al. (2023).

Considering the promising results obtained by local search algorithms recently proposed in the literature (Pérez Cáceres et al., 2019; Cáceres, Araya & Cabrera-Guerrero, 2021; Moyano & Cabrera-Guerrero, 2020; Moyano et al., 2023), in this work, we propose two variable neighbourhood search (VNS) based algorithms to solve the DAO problem. The VNS is a versatile metaheuristic search technique known for effectively exploring large and complex solution spaces, particularly in discrete optimisation problems (Hansen, Mladenović & Pérez, 2010). It is well-suited to this problem due to its structured exploration of multiple neighbourhood structures, which increases the likelihood of escaping local optima and improving solution quality. This systematic exploration of diverse neighbourhoods allows VNS to balance intensification and diversification during the search process (Brimberg, Hansen & Mladenovic, 2010), finding reasonable solutions quickly and making it a robust choice for solving the DAO problem. Specifically, we proposed to use two variations of the VNS algorithms: the variable neighbourhood descent (VND) and the reduced variable neighbourhood search (rVNS). In this study, we evaluate the performance of our proposed algorithms using a set of clinical prostate cancer cases. We benchmark the treatment plans generated by our algorithms against those produced by the traditional sequential method. The findings indicate that our algorithms can produce deliverable treatment plans with a reduced number of apertures and a notable decrease in beam-on time. Additionally, when compared to deliverable plans with an equivalent number of apertures, our approach yields superior results in terms of objective function values. Finally, we compared our algorithm with a local search proposed in the literature (Moyano et al., 2023), obtaining competitive results.

The structure of the article is as follows: In ‘Direct Aperture Optimisation’, we cover the fundamental concepts of IMRT and DAO, along with the mathematical models applied in this research. ‘Variable Neighbourhood Search Algorithms’ describes the algorithms developed and implemented for this study. ‘Computational Experiment’ presents the results from applying our algorithms to a prostate cancer case and discusses the findings. Finally, in ‘Conclusion’, we draw the main conclusions of our work and outline future work.

Direct aperture optimisation

The DAO problem integrates FMO and MLC sequencing, such that the resulting optimisation problem aims to optimise intensities while taking into account the physical limitations of the MLC (Shepard et al., 2002). Although related, this is a completely different optimisation problem as we are no longer interested in optimising beamlet intensities but instead in aperture shapes and their corresponding intensities. This new formulation resulted in a mixed-integer nonlinear model as each beamlet in the apertures becomes a $\left\{0,1\right\}$ (open/closed) variable. The presence of $\left\{0,1\right\}$ variables means the problem is highly complex, and thus, we can no longer rely on complete methods as we did with the FMO problem.

To mathematically represent the IMRT problem, each beam angle is discretised into smaller units called beamlets. At the same time, tissues and the tumour are divided into small sub-volumes known as voxels (Ehrgott et al., 2010). Figure 2 provides a graphical illustration of these concepts.

The IMRT problem can, therefore, be modelled using the representation shown in Fig. 2 (Cabrera-Guerrero et al., 2018a, 2018b). We start by simulating the dosage distribution that is applied to each region’s voxels. The beam angles are separated into $n$ beamlets, where $n$ corresponds to the total number of beamlets for all beam angles, as was previously mentioned. Let $x\in {\mathbb{R}}_{\geqq 0}^n$ be an intensity vector corresponding to $\mathcal{A}$, where $\mathcal{A}$ represents a beam angle configuration (BAC). The vector’s $x_b$ components each indicate how long the $b$-th beamlet was exposed to radiation. Equation (1) (Ehrgott et al., 2010; Cabrera et al., 2018) calculates the radiation dose deposited in voxel $v$ of region $r$ by fluence map $x$.

(1) $$d_v^r(x)=\sum _{b=1}^n(D_{vi}^r{x}_b)\mathrm{\forall }v=1,2,\dots ,m^r.$$

Total voxels in region $r$ is given by $m^r$ in Eq. (1), where $r\in R=\{O_1,\dots ,O_Q,T\}$ indicates an element of the collection of regions. OARs are indexed by $r=O_q$ with $q=1,\dots ,Q$, while the tumour is indexed by $r=T$. $D_{vb}^r\in {\mathbb{R}}^{m^r\times n}$ is the dose deposition matrix for region $r$. The rate at which radiation from beamlet $b$ is deposited into voxel $v$ inside region $r$ is indicated by $D_{vb}^r\geqq 0$ (seen in Fig. 3). When BAC $\mathcal{A}$ is taken into account, the set of all possible intensity vectors is $\mathcal{X}(\mathcal{A}\mathcal{)}\subseteq {\mathbb{R}}^{\mathcal{n}}$.

In this study, we consider Eq. (2) as the objective function of our DAO problem (Romeijn et al., 2003, 2005) Here, parameter $m^r$ is, again, the number of voxels of the region $r$ and $Y_r$ is the desired dose for voxels in the region $r$. The function $(\cdot {)}_{+}$ is the maximum between 0 and $(\cdot )$, $d_v^r(x)$ gives the dose delivered by fluence map $x$ to voxel $v$ of the region $r$ (see Eq. (1)), and ${\underset{\_}{\lambda }}_r$ and ${\overline{\lambda }}_r$ are the penalty weights parameter of under-dose and overdose related to region $r$, respectively. As Eq. (2) shows, we consider a mean-square-error-based objective function in order to obtain actual doses as close as possible to the prescribed ones, penalising those over (under)-irradiated voxels in the OARs (tumour).

(2) $$\underset{x}{min}z(x)=\sum _{r\in R}\left[\frac{1}{m^r}\sum _{i=1}^{m^r}[{\underset{\_}{\lambda }}_r{(Y_r-d_v^r(x))}_{+}^2+{\overline{\lambda }}_r{(d_v^r(x)-Y_r)}_{+}^2]\right].$$

Although the model on the Eq. (2) is convex, we cannot directly apply mathematical programming techniques to obtain the optimal fluence map as we need to account for each aperture shape and its corresponding intensities (obtaining only the optimal intensity vector is not enough to produce a treatment plan that the linear accelerator can deliver). Thus, we need to consider those variables in our model.

Let BAC $\mathcal{A}=\{{\mathcal{A}}_1,\dots ,{\mathcal{A}}_U\}$ be a BAC where the number of beams of such a BAC is $U\in {\mathbb{N}}_{>0}$. Let us now consider a beam angle ${\mathcal{A}}_c$. We then state that the set $\mathbb{H}=\{(P^1,I^1),\dots ,(P^N,I^N)\}$, with $(P^c,I^c)$ being the set of ${\mathrm{\Theta }}^c$ aperture shapes and intensities is a solution to DAO (i.e., a treatment plan that meets deliverable constraints). Each aperture shape $S_i^c\in P^c$ is defined as a binary variable matrix. For a beam angle ${\mathcal{A}}_c$, Fig. 4 provides an example of a tuple $(P^c,I^c)$.

Figure 4 shows the possible values a beamlet can take during the optimisation process. If the beamlet is open (i.e., radiation passes through it), its value equals $1$. If the beamlet is closed (i.e., it is blocked by a MLC leaf), its value equals $0$. Matrix elements equal to $-1$ cover no voxel in the tumour and thus are not considered. It is important to note that the matrix $S_i^c$ is a consecutive $1$’s matrix ( $C1$), i.e., $1$ values in the same row must be consecutive, with no $0$ value in between them. This is due to the MLC physical constraints.

The intensity vector $x$, which is utilised in Eq. (1), must be obtained from the DAO solution in order to evaluate $z(x)$. As a result, we have to add up all of the matrices for every tuple in H. On Eq. (3) a linear combination of the aperture shapes $S_i^c$ and their corresponding intensities $I_i^c$ for beam angle ${\mathcal{A}}_c$ yields the resulting matrix (Moyano et al., 2023):

(3) $$A_c=\sum _{i=1}^{{\mathrm{\Theta }}^c}{S}_i^c\cdot I_i^c.$$

Once we have obtained the aggregated matrix $A_c$, we must map it into an intensity vector $x$. To this end, we correlate the location of each beamlet in the aggregated matrix of beam angle $\mathcal{A}$ to its corresponding position $b$ in the intensity vector $x$ of beam angle $\mathcal{A}$. Figure 5 illustrates the procedure explained above.

Variable neighbourhood search algorithms

VNS algorithms are metaheuristic search techniques that explore different solutions and work exceptionally well in discrete spaces (Hansen, Mladenović & Pérez, 2010). These techniques use different neighbourhood structures to guide the search process and improve the quality of solutions.

As described in the work of Hansen & Mladenovic (2009), the algorithm has different variants. In this article, we decided to use two variations of the VNS algorithms: the VND and the reduced rVNS. These algorithms are explained in ‘Variable Neighbourhood Descent’ and ‘Reduced Variable Neighbourhood Search’, respectively. As in any other Local Search algorithm used on the DAO problem, VNS-based techniques need some key components to be implemented, such as a strategy to generate their initial solution, a method to optimise aperture weights and the neighbourhood definitions. In our case, the proposed VNS-based strategy will use the same components presented in Moyano et al. (2023), which is explained in ‘Initial solution’, ‘Aperature weights optimisation’ and ‘Neighborhood definition’, respectively.

Initial solution

As suggested in Moyano et al. (2023), we employ a set of five predetermined apertures for each beam angle $A_c$ in BAC $\mathcal{A}$ to produce the first solution, which we refer to as the ad-hoc technique. These aperture shapes are shown in Fig. 6. The first aperture has a half-open bottom, the second has a half-open top, the third has a half-open right, the fourth has a half-open left, and the fifth is completely open.

After we have created the ad-hoc (or any other aperture shape), we set the optimal intensity $I_i^c$ for each aperture by using Gurobi solver (Gurobi Optimization, 2020). The model we use to obtain such an optimal intensity vector is presented in ‘Aperture Weight Optimisation’.

Neighbourhood movements

To generate the neighbourhood for both algorithms, we use the two neighbourhood movements proposed by Moyano et al. (2023), called ${\mathcal{N}}_1$, and ${\mathcal{N}}_2$, due to the competitive results obtained by the authors. These movements consider all the scenarios where a leaf can cover a beamlet, so it’s a good approach to getting a good aperture set for DAO. The ${\mathcal{N}}_1$ movement focuses on one single randomly chosen aperture shape. Once this aperture is chosen, for every leaf in it, two neighbours are computed: one that opens the leaf and the other that closes it. Given that, every row is composed of two leaves, i.e., each row has, in general, four neighbours (see Fig. 7). However, as seen in Fig. 8, in some specific cases, we may have only one neighbour per row, depending on the leaf’s current position.

After generating the set of neighbours, we look for possible enhancements in the objective function value. Thus, as proposed in Moyano et al. (2023), we keep moving those leaves that improved our objective function in the same direction until no progress is accomplished. This movement is called the intensification operator (Moyano et al., 2023). We apply this operator until no more improvement is achieved.

Two examples of how the intensification operator works are depicted in Fig. 9 (Moyano et al., 2023). The first example shows how the aperture $b$ was obtained from $a$ by closing the left-hand side leaf to the right. We continue to close the leaf, generating the aperture $c$, as this movement enhance the treatment plan quality. The second example shows how from aperture $d$ we can generate aperture $e$ by opening the bottom right leaf to the right. Since this movement doesn’t lead to an improvement on the treatment’s quality, we end the procedure and go back to aperture $d$.

As noted by Moyano et al. (2023), neighbourhood definitions ${\mathcal{N}}_1$ and ${\mathcal{N}}_2$ result in solutions that are quite similar to the original one. In order to enhance the exploration capabilities of the algorithm, an additional operator, called merge, is included in our LS algorithm. With this additional operator, we aim to provide a new solution that combines all the neighbour solutions that led to an improvement after applying operators ${\mathcal{N}}_1$ and ${\mathcal{N}}_2$. The operation of the merging operator is illustrated in Fig. 10. In this instance, the objective function value was enhanced by three neighbouring solutions. After that, we combine these three apertures to create a single new aperture that incorporates the movements produced at each neighbour independently. For a more detailed explanation on the merge operator, please see Moyano et al. (2023).

The main difference between movements ${\mathcal{N}}_1$ and ${\mathcal{N}}_2$ is that the leaves that are modified (either by opening or closing them) belongs to different aperture shapes (movement ${\mathcal{N}}_1$ is applied to the same randomly chosen aperture shape). Just as in movement ${\mathcal{N}}_1$, we create one or two neighbours for every leaf, just like in the ${\mathcal{N}}_1$, as seen in Fig. 11. We also apply the intensification and the merge operators to the obtained solutions after applying ${\mathcal{N}}_2$.

Computational experiment

This section introduces the experiments performed using the VND and the rVNS algorithms and analyses the results obtained. ‘Experimental Setup’ describes the prostate case used in the experiment. Then, in ‘Experiments’, the analysis of the algorithm is made. To leverage the two neighbourhood moves ( ${\mathcal{N}}_1$, and ${\mathcal{N}}_2$) in our VND algorithm, we have defined two variants: $VN{D}_{12}$ and $VN{D}_{21}$. These variants differ in the order in which the neighbourhood moves are applied. Specifically, $VN{D}_{12}$ applies move ${\mathcal{N}}_1$ first, followed by move ${\mathcal{N}}_2$, meanwhile $VN{D}_{21}$ applies move ${\mathcal{N}}_2$ first, followed by move ${\mathcal{N}}_1$. Finally, in ‘Clinical Indicators Analysis’, we present a clinical indicator analysis comparing the quality of the treatment plans generated by our rVNS algorithm with those obtained using the DAO implementation available in the matRad platform (Wieser et al., 2017).

Experimental setup

This study uses a prostate case from CERR package (Deasy, Blanco & Clark, 2003) to evaluate our algorithms’ performance. In addition, we used four cases (TRT001, TRT002, TRT004 and TRT005) of the CAS dataset (Cabrera-Guerrero et al., 2018a) to evaluate the performance of the clinical criteria of our best algorithm with the method proposed by matRad. These cases consider two organs at risk, namely the bladder and the rectum, as well as the prostate (PTV) (see Fig. 12).

For the CERR case study, the prostate has 15,172 voxels, the bladder has 22,936 voxels, and the rectum has 18,128 voxels. The cases of CAS, as shown in Table 1, range from 9,806 to 19,652 voxels for the tumour, 19,762 to 52,024 voxels for the bladder, and 8,500 to 29,346 voxels for the rectum. Regarding computational cost and scalability, it is important to note that while the rVNS framework itself demonstrates good scalability to more anatomically complex cases, solver performance may be impaired as the number of regions of interest increases. However, complex anatomical sites, such as head-and-neck, typically present smaller target volumes and OARs compared to prostate cases, which can reduce the overall problem size and mitigate this effect.

Table 1:

Voxels per region in CAS cases.

Case	Prostate	Bladder	Rectum
TRT001	13,081	19,762	8,500
TRT002	19,652	52,024	10,221
TRT004	9,806	22,597	12,849
TRT005	10,082	23,684	21,171

DOI: 10.7717/peerj-cs.3094/table-1

We consider 72 available beam angles, all on the same plane in all the cases. Just as in Moyano et al. (2023), we consider a set of 14 equidistant BACs to do our experiments, as is shown in Table 2. We run our algorithm 30 times per BAC, being 30 a widely accepted value for statistical analysis (Hays & Winkler, 1970). On average, the methods need around 10 min and depend mainly on the size of the prostate case (number of voxels and beamlets) as well as on the quality of the initial set of apertures. Thus, we acknowledge that better times can be achieved if more effective initial apertures are generated. This is part of our future work. The maximum memory required by the algorithms was 6.5 GB.

Table 2:

Initial beam angle configurations.

BAC	${\mathit{\theta }}_{\mathbf{1}}$	${\mathit{\theta }}_{\mathbf{2}}$	${\mathit{\theta }}_{\mathbf{3}}$	${\mathit{\theta }}_{\mathbf{4}}$	${\mathit{\theta }}_{\mathbf{5}}$
1	0	70	140	210	280
2	5	75	145	215	285
3	10	80	150	220	290
4	15	85	155	225	295
5	20	90	160	230	300
6	25	95	165	235	305
7	30	100	170	240	310
8	35	105	175	245	315
9	40	110	180	250	320
10	45	115	185	255	325
11	50	120	190	260	330
12	55	125	195	265	335
13	60	130	200	270	340
14	65	135	205	275	345

DOI: 10.7717/peerj-cs.3094/table-2

The prescribed dose values and associated penalty weights used in the objective function $z(x)$ (see Eq. (2)) are presented in Tables 3 and 4. Table 3 corresponds to the same parameters used in Moyano et al. (2023), enabling direct comparison of results. Table 4, on the other hand, presents the parameters used for the CAS dataset, where more stringent dose limits were imposed on OARs to reflect a more clinically demanding scenario.

Table 3:

Value of ${\mathit{T}}_{\mathit{i}}$, ${\underset{―}{\lambda }}_i$ and $\overline{{\mathit{\lambda }}_i}$ for function z(x) CERR case.

Organ	${\mathit{Y}}_{\mathit{r}}$	${\underset{―}{\lambda }}_r$	$\overline{{\mathit{\lambda }}_r}$
PTV	76 Gy	5	5
Rectum	65 Gy	0	1
Bladder	65 Gy	0	1

DOI: 10.7717/peerj-cs.3094/table-3

Table 4:

Value of ${\mathit{T}}_{\mathit{i}}$, ${\underset{―}{\lambda }}_i$ and $\overline{{\mathit{\lambda }}_i}$ for function z(x) for CAS cases.

Organ	${\mathit{Y}}_{\mathit{r}}$	${\underset{―}{\lambda }}_r$	$\overline{{\mathit{\lambda }}_r}$
PTV	76 Gy	5	5
Rectum	55 Gy	0	1
Bladder	50 Gy	0	1

DOI: 10.7717/peerj-cs.3094/table-4

Tables 3 and 4 shows the penalty weights for under/over irradiated voxels, ${\underset{\_}{\lambda }}_r$ and ${\overline{\lambda }}_r$, respectively. Also, we consider an intensity range of $\left[0,20\right]$ for ${\gamma }^{min}$ and ${\gamma }^{max}$, respectively. Finally, we consider five apertures for each beam angle $A^c$ for both datasets.

Conclusion

This article compares different algorithms based on variable neighbourhood search combined with mathematical programming to generate a radiation therapy treatment plan for cancer by solving the DAO problem. The methods proposed in this article can identify a range of aperture shapes and intensities for each beam angle of a specified BAC. Unlike the ones obtained by the traditional sequential approach, all these treatment plans are deliverable, and they all exhibit clinically acceptable delivery times. Furthermore, despite the algorithm’s limitation to five aperture shapes per beam angle, resulting in a maximum of twenty-five aperture shapes, our proposed methods can generate treatment plans that yield highly competitive objective function values.

The algorithms discussed in this article achieved reduced treatment planning time compared to traditional sequential methods, allowing a quicker plan generation without sacrificing quality. This is crucial for enhancing patient outcomes, as better-optimised plans can reduce side effects and improve the accuracy of tumour targeting, ultimately contributing to more effective cancer control. The ability to generate competitive treatment plans even under constraints, such as using only a limited number of apertures per beam, highlights the practicality and applicability of these algorithms in clinical settings, potentially leading to improvements in both treatment effectiveness and patient safety.

While the results obtained in this study demonstrate the effectiveness of the rVNS proposed algorithm compared with the LS method reported in the literature, it is important to highlight certain limitations. Although the rVNS algorithm does not outperform the matRad DAO tool, it still produces competitive treatment plans from a clinical indicators perspective. These results motivate further improvements to our method, particularly in balancing the trade-off between target coverage and OAR sparing. Importantly, our current implementation restricts the solution to only five aperture shapes at most, which limits the algorithm’s flexibility. As highlighted in Cáceres, Araya & Cabrera-Guerrero (2021), reducing the number of apertures can significantly degrade plan quality. Since matRad has no such constraint, it can more effectively conform to the PTV while sparing the OARs, especially in anatomically complex scenarios. However, this better coverage achieved by matRad comes at a cost of more aperture shapes, which in turn lead to more MU delivered to patients. Thus, while matRad favours the PTV coverage at a cost of more complex treatment plans, our algorithm obtains quite competitive results, keeping the complexity of the treatment plan low and limiting the number of MU of the produced treatment plans.

While this study focused on prostate cancer cases, which provide a well-established starting point for DAO research due to their regular geometry and common use in benchmark datasets, we acknowledge that validating our method on more anatomically complex tumour sites, such as head-and-neck cases, is essential to assess its broader applicability. Our local search strategy, based on small, incremental neighbourhood movements, is inherently adaptable and should generalise well to irregular anatomical configurations, though.

In addition, we can identify various research directions to improve the results obtained by our approach. First, exploring new neighbourhood movements would allow us to find better-quality treatment plans. Some interesting movements that can be adapted from other algorithms used in the literature include the crossover movement of the genetic algorithm proposed by Cao et al. (2009) and the repair heuristic used for the particle swarm optimisation in Tello-Valenzuela, Moyano & Cabrera-Guerrero (2023). Furthermore, other techniques, such as memetic algorithms, can be explored to leverage the ability of LS to find high-quality aperture shapes and the genetic algorithms’ capacity to explore the search space. Finally, focusing on new strategies to generate better initial solutions within our framework might be a worthwhile research direction, as this would enable the algorithm to further reduce computational time.

Finally, regarding computational cost and scalability, it is worth noting that the rVNS framework exhibits good adaptability to more anatomically complex cases. While solver performance may be affected as the number of regions of interest increases, sites such as head and neck typically involve smaller target and OAR volumes, which can mitigate this effect. Moreover, the tendency of our algorithm to favour solutions with fewer active apertures further contributes to maintaining efficient computation times. As part of our future work, we plan to conduct a more detailed analysis of runtime behaviour across a broader range of anatomical scenarios to better support the clinical applicability of the proposed methods.