Abstract
One important problem in radiation therapy for cancer treatment is the selection of the set of beam angles radiation will be delivered from. A primary goal in this problem is to find a beam angle configuration (BAC) that leads to a clinically acceptable treatment plan. Further, this process must be done within clinically acceptable times. Since the problem of selecting beam angles in radiation therapy is known to be extremely hard to solve as well as time-consuming, both exact algorithms and population-based heuristics might not be suitable to solve this problem. In this paper, we compare two matheuristic methods based on local search algorithms, to approximately solve the beam angle optimisation problem (BAO). Although the steepest descent algorithm is able to find locally optimal BACs for the BAO problem, it takes too long before convergence, which is not acceptable in clinical practice. Thus, we propose to use a next descent algorithm that converges quickly to good quality solutions although no (local) optimality guarantee is given. We apply our two matheuristic methods on a prostate case which considers two organs at risk, namely, the rectum and the bladder. Results show that the matheuristic algorithm based on the next descent local search is able to quickly find solutions as good as the ones found by the steepest descent algorithm.
1. Introduction
Radiation is one of the most common therapies used to treat patients suffering from cancer. The purpose of radiation therapy is to deliver a dose of radiation to a tumour in order to sterilize all cancer cells and to minimize the collateral effects on the surrounding healthy organs and tissues. Intensity modulated radiation therapy (IMRT) is the most common technique within radiation therapy.
We can separate IMRT problem into three sequential optimisation (sub-)problems: the beam angle optimisation (BAO) problem, the fluence map optimisation (FMO) problem, and the multileaf collimator sequencing problem, Ehrgott et al. [1]. In the BAO problem, we determine the number and directions of the beam angles we shall use to produce a treatment plan. The set of beams used to treat a patient is called beam angle configuration (BAC). Then, in the FMO problem, we determine the radiation intensities that will be delivered from each beam angle. The solution to this problem is a vector of intensities that is called fluence map. Finally, a sequence of movements of a physical device called multileaf collimator is computed in order to efficiently deliver the fluence map computed during the previous phases. It is clear from here that the selection of beam angles in the BAO phase has a big impact on the quality of the fluence map computed in the FMO phase. That is, a good combination of beam angles will lead to a good quality fluence map and, consequently, will produce a good quality treatment plan. In this paper the problem of selecting a good quality BAC is addressed.
To measure the quality of a BAC we need to solve the associated FMO, that is, we solve the FMO problem for each evaluated BAC. Computing the optimal fluence map for a BAC is a time-consuming process. Further, one practical constraint when solving the BAO problem is that we need to obtain results within clinically acceptable times. In our experience, “clinically acceptable times” are around 12 hours (i.e., running algorithms overnight) so treatment planners can decide among alternative treatment plans during the day after patients images have been obtained. Thus, we need to find efficient strategies that produce good quality treatment plans within these time limits; that is, not too many BACs can be evaluated during the optimisation process. For this reason, sophisticated (meta-)heuristic algorithms such as population-based algorithms might not be suitable for solving this problem.
The aim of this paper is to compare two different local search strategies to capture the trade-off between solutions quality and the required time before convergence. Using the next descent algorithm, we aim to accelerate the search without any major impairment to treatment plans’ quality.
The remainder of this paper is organised as follows. In Section 2 we introduce the BAO problem we aim to solve in this paper and a brief review on different approaches dealing with the BAO problem is presented. In Section 3 both the steepest descent and the next descent algorithms are presented and their main differences are highlighted. In Section 4 the set of instances considered in our experiments are presented. Obtained results are also discussed in this section. Finally, in Section 5 some conclusions are drawn and the future work is outlined.
2. Intensity Modulated Radiation Therapy: An Overview
During the last three decades or so, many researchers have worked on the problem of finding efficient treatment plans for radiation therapy for cancer treatment. Most of their efforts have been focused on the problem of finding the optimal fluence map that can be delivered to the patients given a predefined BAC (see, e.g., [2–5]). Unfortunately, significantly less attention has been paid to the problem of finding the best BAC. This is, in part, because the BAO problem is, from a mathematical point of view, much harder to solve than the FMO problem. Further, treatment planners rely on their experience to choose BACs that are good enough to produce clinically acceptable treatment plans. Figure 1 shows the traditional try-and-error process followed by treatment planners in clinical practice.
As we can see, treatment planners try an initial BAC based on their experience and judgement. If no clinically acceptable treatment plan can be produced using the initial BAC, then the expert proposes a new BAC to be tested. This process is repeated until a clinically acceptable treatment plan is produced. Depending on the particular case and the treatment planner experience, this process can take very long.
The first attempt to study the BAO problem from a theoretical point of view is presented in Bortfeld and Schlegel [6]. Since then, both exact and heuristic methods have been proposed to solve this problem. For instance, Rocha et al. [7] propose a guided pattern search method to solve the BAO problem. Genetic algorithms have also been used to explore the beam angles space in Dias et al. [8]; Lei and Li [9]; and Li et al. [10]. In Li et al. [11], the authors propose an ant colony optimisation algorithm that uses a similar representation to Li et al. [10]. In Li and Yao [3], authors propose a hybrid algorithm that combines ant colony optimisation and genetic algorithms. The authors claim that their approach is faster than previously implemented genetic algorithms. A particle swarm optimisation algorithm has also been proposed in Li et al. [12]. In Aleman et al. [13], the authors propose a surface response to find high-quality BACs. We need to highlight at this point that the population-based algorithms as well as the evolutionary algorithms cited above address the BAO problem with only linear or quadratic objective functions and constraints. Thus, they can evaluate a large number of BACs within time limits. This is not the case for our study, as the objective function and constraints we consider here are highly complex nonlinear functions. We will explain this at the end of this section.
Local search approaches have also been proposed. For instance, in Aleman et al. [14] local search algorithms, namely, simulated annealing and add/drop, are implemented. Authors highlight the fact that good quality BACs can be generated using this kind of strategies. Further, they also propose a neighbourhood structure that allows for faster convergence, reducing the time required to obtain a good quality plan. Simulated annealing is also used in Mohammadi et al. [15] and Dias et al. [16].
Methods that combine heuristic and mathematical programming have also been proposed. For instance, a novel method that combines simulated annealing with gradient-based search is presented in Bertsimas et al. [17]. Similarly, in Craft [18], authors combine a local search algorithm with linear programming and gradient search. In Lin et al. [19], authors propose a two-stage heuristic method where, in the first stage, Benders decomposition is considered. Obtained solutions are then deputed during the second stage using local branching.
In Cabrera et al. [20] the authors propose a two-phase approach to solve the multiobjective version of the BAO problem. One distinctive feature of their approach is that they use a single objective local search during the first phase to produce efficient treatment plans.
In all the approaches above, the quality of a BAC is determined by the best fluence map such a BAC is able to produce. Thus, for each BAC we want to evaluate, we must solve its associated FMO problem. To mathematically model the BAO problem, each beam angle is discretised into a set of beamlets that radiation is delivered from. Similarly, organs at risk (OARs) and the tumour are discretised into a set of voxels. Figure 2 shows a representation of the problem considering two beam angles (labelled as beam), one tumour (labelled as Tumour), and two organs at risk (labelled as Normal Tissue). As we can see, each beam is discretised forming a “grid” of beamlets (in blue), while the tumour and OARs are divided into voxels (in orange and green, resp.)
Although all beam angles around the patient body are pointing to the centre of the tumour, we cannot avoid the fact that the radiation delivered from beamlets hits not only the tumour but also some OARs (see Figure 3). Thus, finding fluence maps that generate uniform doses on the tumour and that avoid as much as possible overirradiating the OARs is a key task during the treatment plan generation process.
A fluence map specifies the total time each beamlet is allowed to irradiate the tumour. Figure 4 shows an example of a fluence map for a beam angle. In this illustration, the beam angle is divided into a grid of beamlets. Using a multileaf collimator (a physical device consisting of a set of leaves that can block the radiation coming from the linear accelerator), radiation can be modulated so we can determine how much radiation is delivered from each beamlet. Radiation delivered from each beamlet is represented as a bar in Figure 4, where higher bars mean more radiation. While some beamlets will deliver no radiation during the treatment plan, others will irradiate not only the tumour but also some OARs. Discretizing beam angles into beamlets and organs into voxels allows modelling the problem as an optimisation problem. We first let be the set of beam angles that radiation can be delivered from. Although in this paper we only consider coplanar angles, our algorithms can be easily extended to noncoplanar angles. Angular resolutions other than can also be used. Let be a BAC where is the set of all -element subsets of , with being the number of beam angles we shall consider within a BAC. We denote the th angle of BAC by , for .
As mentioned before, in the BAO problem we need to solve, for each evaluated BAC, the associated FMO problem. Thus, we first introduce the mathematical model of the FMO problem we consider in this work for a fixed BAC :with , where is an objective function and the fluence map is a nonnegative vector with each component , called fluence, representing the length of time that a patient is exposed to beamlet and where is the total number of beamlets summed over all possible beam angles. The set is the set of all feasible solutions of the FMO problem when BAC is considered. Only beamlets that belong to a beam angle in are allowed to have a value greater than zero.
Different objective functions have been proposed in the literature to solve the FMO problem. The vast majority of this functions need to compute the total radiation dose deposited into each voxel of the tumour and OARs by fluence map . This dose delivered on region by a fluence map is represented by a vector we call , where is the total number of voxels in region and each of its elements is calculated as in (2) (see [1]).where is an element of the index set of regions where refers to the tumour and refers to the -th OAR with . Matrix , called dose deposition matrix, is a given matrix where its elements determine how much radiation is delivered for each beamlet to each voxel per time unit at region .
In this paper we consider an objective function originally proposed by Wu et al. [4] and called the logistic function. This function is based on a well-known measure used in radiation therapy called the generalised equivalent uniform dose (gEUD) which was proposed in the work by Niemierko [21]. The gEUD is defined in the work by Niemierko [21] as “the equivalent dose, which if distributed uniformly across the tumour, would lead to the same level of cell killing as the actual dose distribution of interest.” Mathematically, the gEUD is as follows:where is a region dependant parameter defined by treatment planners.
The logistic function originally proposed in Wu et al. [4] to compute the optimal fluence map for an associated BAC iswhere
Parameters and are the prescribed gEUD value for tumour and the maximum gEUD value allowed for OARs, respectively, and are user-defined parameters that indicate the importance of the tumour and the -th OAR, respectively. gEUD-based objective functions as the ones in (4) have been largely used to solve the FMO problem (e.g., [4, 5, 20, 22]).
In this paper we slightly modify the logistic function in (4) by removing (5) from the objective function and adding a constraint to the model setting gEUD of the tumour to , as proposed in Cabrera et al. [20]. Thus, objective function of the FMO problem we consider in this paper is as follows:As discussed in the work by Cabrera et al. [20], the is a convex function of and is also convex in for the region of interest. Thus, problem (7) is convex and, therefore, optimal solutions can be found by exact algorithms that can deal with nonlinear problems.
Having defined the FMO problem, we can now introduce the mathematical model for the BAO problem we aim to solve in this study.
The BAO problem has been shown to be nonconvex and highly nonlinear. As an example of this, Figure 5 shows a surface plot for the logistic function in (7) for all possible BACs that consider two beam angles. As we can see, the surface presents several combinations of beams where locally optimal solutions can be found. Clearly, as more beam angles are added, the problem becomes much harder to solve. Also, we can note that worst values are obtained when both angles within the BAC are the same (or too close to each other). Opposite angles, that is, with a difference around 180°, also obtain high values. We will discuss this point in the conclusions of this study.
3. Local Search-Based Matheuristic Algorithms for BAO
As mentioned before, many heuristic algorithms have been proposed in the literature to solve the BAO problem. In this paper we propose to solve this problem by using local search algorithms. First, we implement a simple yet efficient local search algorithm called steepest descent, which was previously used in Cabrera et al. [20] to solve the multiobjective version of the BAO problem. Algorithm 1 shows the pseudocode of this method.
| ||||||||||||||||||||||||||||
The steepest descent algorithm starts with an initial BAC, which can be either provided by the treatment planner or randomly generated. This initial solution is labelled as the current solution. For this current solution, a set of neighbours is generated by applying the neighbourhood move (as explained in the paragraphs below) on the current solution. Once the entire neighbourhood defined by has been generated, the best neighbourhood with respect to its objective function is selected. If the best neighbour solution is better than the current one, then the best neighbour is labelled as the new current solution. If the best neighbour is not better than the current solution the algorithm stops and the current solution is returned as the locally optimal solution the algorithm found. As before, the steepest descent algorithm is very simple to implement and it converges to good quality solutions. However, since it evaluates the entire neighbourhood of the current solution at each iteration, sometimes it can take too long to converge. Since one mandatory requirement in clinical practice in radiation therapy is to provide good quality treatment plans within clinically acceptable times, algorithms such as the steepest descent might not be suitable. Thus, in this paper, we propose an algorithm that is able to produce treatment plans faster than the steepest descent without any major impact on the quality of the obtained treatment plans. Algorithm 2 shows this algorithm called the next descent algorithm. Just as in the steepest descent algorithm, the next descent algorithm starts with an initial BAC which is also labelled as the current solution. Then, for each element within the neighbourhood of the current solution the FMO problem is solved and the quality of the resulting treatment plan is compared to the current solution. If the neighbour is better than the current solution, then the neighbour is labelled as the new current solution and its neighbourhood is generated.
| ||||||||||||||||||||||||||||||||
If the current solution is better than the neighbour solution, then the neighbour solution is removed from the neighbourhood and the next neighbour is evaluated. The algorithm stops once the neighbourhood is empty; that is, no neighbour solution is better than the current solution. Since the next descent algorithm does not explore the entire neighbourhood of the current solution at each iteration, it needs fewer objective function evaluations, leading to a faster convergence.
Both the steepest descent and the next descent algorithms need a neighbourhood to be defined. In this paper we consider the following neighbourhood:where, again, . Given this neighbourhood definition, we know that neighbours will be generated. Further, in some cases, we can have that with . We allow neighbours to repeat a beam angle although it is well known that this situation impacts the neighbour quality. We need to point out that other neighbourhood definitions that support noncoplanar beam angles can be found in Lim and Cao [23]; Mišić et al. [24]; Yarmand and Craft [25].
4. Computational Experiments
We present in this section the computational experiments carried out in this work. We use an Intel i7 processor and 32 GB of RAM to run our experiments. Linux 16.04 was the operating system. Both the steepest descent and the next descent algorithms were coded in JAVA 8 language using NetBeans IDE. Experiments are performed on a prostate obtained from the Computational Environment for Radiotherapy Research (CERR, http://www.cerr.info/). The tumour is located at the prostate while OARs considered in this study are the bladder and the rectum. We need to point out that prostate cases are, in general, very difficult as the tumour and OARs are, as in our case, overlapped (see Figure 6).
In this prostate case, the tumour has more than 7,000 voxels, the rectum about 5,500, and the bladder around 9,500. Although the same multileaf collimator is used for all beam angles the number of beamlets considered in the optimisation process depends on each beam angle. This is because those beamlets that do not irradiate at least one voxel of the tumour are not considered as decision variables in the optimisation process. Thus, the number of decision variables (beamlets) depends on the BAC and ranges between 160 and 220. The number of beam angles considered in a BAC is equal to 5. As we mentioned before, we use IPOPT Wächter and Biegler [26] as the solver for all nonlinear optimisation problems.
Table 1 shows the value of the parameters of the logistic function used for our experiments. These values are the same as in Cabrera et al. [20].
Three sets of initial BACs are generated for our experiments. The first set () is a set of 14 BACS for which beam angles are equidistant. The second set () is a set of initial BACs for which beam angles have been selected randomly but considering some geometrical constraints such as avoiding opposite beam angles and angles that are too close to each other. The third set of initial BACs () consists of BACs for which beam angles have been selected absolutely randomly. All BACs evaluated in our experiments consist of 5 beam angles. We need to point out that for the SD algorithm only one run per initial BAC is needed, as it is a deterministic algorithm. For the ND, 10 independent runs are performed for each initial BAC. Reported results in the next section correspond to the average of these 10 runs for each initial BAC.
4.1. Results
A summary of the obtained results is presented in this section. Tables 2, 3, and 4 show the results obtained by both the steepest descent (SD) and the next descent (ND) algorithms for instance sets , , and , respectively. As we can see, not only is the next descent algorithm faster than the steepest descent algorithm in all but three instances in total but also it performs slightly better, in average, than the steepest descent algorithm in sets and . For set , the steepest descent performs slightly better than the next descent algorithm. These results confirm that we can use faster algorithms without any major impairment in the treatment plan quality.
Further, using faster algorithms such as the next descent algorithm will allow us to, for instance, explore different parts of the beam angle space by restarting the algorithm from different initial BACs. Moreover, this type of algorithms might also be used within more complex frameworks that combine different optimisation algorithms.
Figure 7 shows the % of improvement (savings) of the next descent algorithm compared to the steepest descent algorithm for the all 44 instances evaluated in this study. Positive values mean next descent performs better. Negative values mean steepest descent performs better.
Similar to Figure 7, Figure 8 shows the variations in the number of iterations each algorithm needs to perform before convergence. Again, positive values mean that next descent algorithm was faster. It is interesting to note that, for some few instances (, , and ), the next descent algorithm took longer (4,897; 1,280; and 6,483 seconds, resp.) than the steepest descent. Although possible, this is very unusual.
As a summary, results show that the next descent algorithm we implement in this paper is able to obtain treatment plans that are quite competitive with respect to the ones obtained by the steepest descent algorithm. The next descent algorithm takes a fraction of the time that the steepest descent takes, though.
5. Conclusions
In this paper we have implemented two local search methods, namely, the steepest descent and the next descent algorithms. While the former needs to explore the entire neighbourhood at each iteration, the latter does not. As a consequence, the next descent algorithm converges much faster than the next descent algorithm, as it explores a smaller part of the solution space.
We try these two local search methods on a complex problem arising in radiation therapy for cancer treatment called the beam angle optimisation problem. Both algorithms are able to produce clinically acceptable treatment plans. However, the next descent algorithm converges to locally optimal solutions much faster than the steepest descent. Results show that the quality of the treatment plans obtained by the steepest descent is, in average, slightly better than the ones obtained by the next descent algorithm. However, for almost half of the instances we test in this paper, the next descent algorithm demonstrates to perform better than the steepest descent.
As a future work, other local search algorithms such as Tabu Search and Variable Neighbourhood Search must be tried with the aim of obtaining high-quality treatment plans within clinically acceptable times. Moreover, the multiobjective version of the next descent algorithm can be tested on the MO-BAO problem.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
Guillermo Cabrera-Guerrero wishes to acknowledge both DI 039.416/2017 Project and CONICYT/FONDECYT/INICIACION/11170456 Project for partially supporting this research.