A matheuristic for the selection of beam directions and dose distribution in Radiotherapy Planning

In this paper a matheuristic using a combined Genetic Algorithm (GA) and exact method approach is proposed for selecting the position of the beams and dose distribution in Intensity Modulated Radiotherapy Planning (IMRT). GA selects a set of beams, for which the dose distribution is determined in the process of the GA's evaluation, using an optimisation model that is solved by an Interior Point method. Two instances are used to evaluate the performance of the matheuristic, comparing to the optimum solution, in terms of solution and computation time, found using the exact methodology of Branch and Bound. The results show that the matheuristic is appropriate to this problem in the case study proposed, as it is extremely faster than the exact method and also have reached the optimum solution in several experiments done.


Introduction
Radiotherapy is one of the most used kinds of treatment to combat cancer. The idea of the treatment is to determine the necessary dose to reduce, or even kill the tumor, whilst making sure not to affect the surrounding healthy organs.
One of the first optimisation models for the radiotherapy problem was proposed by Bahr et al. [3], in 1968. Thereafter, many improvements have taken place in the fields of radiotherapy and Operational Research, to improve the treatment planning.
However, the success of the treatment planning depends on the experience of the physicist in using a trial-error procedure to solve the following problems.
• How to position the gantry and consequently the beam direction (geometric problem); • How much dosage to deliver through each beam and sub-beam (dose intensity problem); These are classical problems in radiotherapy planning and the literature reports some papers that use either exact methodologies or metaheuristics to solve them.
One of the most employed exact method to solve these problems is Interior Point. It was used by Obal et al. [19] to solve the multiobjective model of the dose intensity problem. Romeijn et al. [21] proposed a linear approach for the same problem and discussed the hard restrictions of the problem.
Viana et al. [27] used Interior Point to solve the model proposed by [1] -to the dose intensity problem -adding an heterogeneity correction in the composition of different types of irradiated tissues, based on proportions among their different linear attenuation coefficients.
Voet [32] and Breedveld [5] used the algorithm iCycle, that is an algorithm for integrated, multi-criteria optimization of beam angles and profiles, based on a priori defined plan criteria, that is called a "wish-list", and is solved using Interior Point.
Other approaches use Branch and Bound or Branch and Cut methodologies [17], [25] and [12]. Souza [25] developed an integer optimisation model to solve the beam directions problem and to evaluate the average dose in the organs at risk. The model was solved using a Branch and Cut algorithm. Lim et al. [17] developed a computational approach to optimise the beam directions and fluence map in IMRT. The authors use an iterative procedure trying to improve the computational performance in comparison with the integer model that was solved with the Branch and Cut methodology.
Gevert et al. [12] proposed a single model for the optimisation of the beam directions and the dose intensity. A set of candidate beams is proposed, and the model inserts an upper bound on the number of beams that could be chosen. The choice is made to optimise the delivery dose in each tissue, finding the distribution of the dose in the beams.
The approach to select a limited number of beam directions is very important, because it determines the gantry position, taking into consideration the limited number of beams imposed by the health plan. That is a clinic reality in the Brazilian health system, although is not in some health institutions in other countries, as discussed in [6].
Moreover, some research has used metaheuristics for radiotherapy problems. Petrovic et al. [20] used a case based reasoning system (CBR), that employs a modified DempsterShafer theory to fuse dose plans suggested by the most similar cases retrieved from the case base, intending to find the best selection of dose distribution and beam directions. After each run of the CBR system, the Simulated Annealing algorithm searches for such a combination of weights for which the difference between the doses suggested by the oncologist/system and doses recommended by the extracted similar cases are minimal.
Simulated Anealling is used by Acosta et al. [1], who propose a model that aims to determine the optimum distribution to each radiation beam, and by Aubry et al. [2].
Another approach was proposed by Hou et al. [15] who developed an optimum solution process in two phases: (1) the optimization of the profiles to a certain set of beams; and (2) the selection of an ideal set of beams based on the result of the objective function to the optimization of the profile. A Dynamic Simulation Algorithm was used to solve the first phase and Genetic Algorithms was used in the second phase.
The research proposed in this paper is a matheuristic to solve a goal programming, which is a hybridization between Genetic Algorithm (GA) and Interior Point. GA is proposed to select the beam directions, and Interior Point to solve the dose distribution problem to GA's selected set of beams.
The experiments have been conducted to improve the computation time that the methodology of Gevert et al. [12]. Although their methodology is robust, it could be very expensive in computation time, which is not desirable due the urgent natural of radiotherapy planning.

Problem formulation
Conventionally, the IMRT planning is done on a linear accelerator, with the following components: the couch, the gantry and the collimator (Figure 1).
The linear accelerator is an isocentral machine, that means there is a system formed to the plane of the gantry (that could rotated 360 degrees) and the plane of the couch. The central point of this system is called the mechanical isocentre ( Figure 1). This point is first marked on the computational tomography (CT) slices of the patient, and after that, some tags are placed on the body of the patient that serve as reference to position the couch and the gantry. The mechanical isocentre is the base point that represents the centre of the sphere where it is possible to position the beam direction. The problem is modeled considering a region of the body obtained by the slice of the CT, as shown in Figure 2. The structures selected are: the tumor, in the prostate region; the organs at risk (noble tissue), represented by the head of femur, rectum and the bladder. The CT slice is represented by a grid with p pixels, in which each pixel is considered as healthy, noble or tumor tissue.
The central objective of radiotherapy planning is to determine how to design the fluence map so that the dose that reaches the healthy and noble pixels is as small as possible and that the tumor dose is closest to the one prescribed by the Physician. This is a goal programming problem [16], because it has two clear set of goals: delivery of the necessary dose to the tumor, and least dose possible to the healthy organs. These are mathematical goals and deviation variables should be used, as will be shown in section 3.
The objective of this paper is to develop a matheuristic to find the beam directions and the dose intensity, intending to achieve the treatment goals.

Model
Attending to the treatment goals, model (1) is proposed in Obal et al. [19]. This model aims to find the best distribution of the dose for an established set of beams.
pn, ph and pt are the respective number of noble, healthy and tumor pixels, with pn+ph+pt = l.c, and l × c is the CT slice's size. X ∈ R K.m.n represents the emitted dose by each sub-beam, in which m.n is the number of sub-beams in the beam k, with k = {1, ..., K}.
The dose delivered to the noble and healthy tissue needs to respect the dose limit: S n ∈ R and S s ∈ R, respectively. D ∈ R is the dose intensity prescribed by the doctor, that needs to be in the tumor structure. B n ∈ R pn , B s ∈ R ph and B t ∈ R pt , index each kind of tissue: noble, healthy and tumor, respectively.
A n ∈ R pn×m.n , A s ∈ R ph×m.n and A t ∈ R pt×m.n represent the absorption dose in each tissue, by the emitted dose of each sub-beam.
The vectors: θ + ∈ R pn , δ + ∈ R ph and + ∈ R pt represent the deviation of excess dose in the noble, healthy and tumor tissues, respectively; − ∈ R pt represent the deviation of deficit dose on the tumor structure; θ − ∈ R pn and δ − ∈ R ph are relate to the difference between the dose that arrives in the noble and healthy tissues and its respective bounds. So, this model has 2.(pn + ph + pt) + K.m.n variables and l.c constraints, and aims to find as closer as possible dose D in the tumor structures, which is represented on the functions (c) and (d), and the minimum dose in the noble and health structures, represented in the functions (a) and (b).
The model (1) include the constrains (2) and (3), which aims to determine which beams are optimum to be selected and how much intensity dose in each selected beam is best achieve the optimal radiotherapy plan.
in which K is the possible number of beam direction; η is the maximum number to be select, with η ≤ K; z k is a binary variable that represent if the k beam is used or not; and X k ∈ R m.n is the sub-vector of X, that represent just the beam k. The addition of constraints (2) and (3) to the model (1), results in a mixed integer linear programming model. This leads to more computation time of the exact solution process, because of the problem dimension, as will be show in section 5.1. So that, the Genetic Algorithm (GA) is proposed in this research, that make the selection of the beams to be used.

Matheuristic for selection of beam directions and dose distribution
This research proposes a hybrid methodology, that uses GA and mathematical programming exact method. For each set of beams that GA proposes, the dose distribution is determine by an Interior Point method solution to model (1).
Based on the principles of natural selection and genetics [11], [4] [14], the Genetic Algorithm is a metaheuristic that generates near-optimal solutions for optimisation problems. It is used due to its ability to find pragmatic solutions in reasonable computation time.
In this algorithm we propose the individual structure as an array of dimension 1 × K where K is the proposed number of possible angular measurements or beam positions. Each column (gene) carries the information 0 or 1 that represent if that beam direction will be used or not, respectively. Figure 3 illustrates the structure of an individual with 8 possible angles at which the first, third, sixth, seventh and eighth beams were selected to make up the treatment plan. The initial population is generated using a specific number of individuals. All the individuals created in the initialization are feasible and represent the set of beam directions used in the planning.
Then, this initial population is improved using the following steps [10]: (i) Fitness evaluation: The best set of beams is the set which can provide the best possible treatment. The treatment quality is a combination of the set of beams and the dose distribution that it can provide. These problems are dependent, because a better dose distribution depends on the beams selected. As the beam selection is determined by the GA's individual, the dose distribution problem is solved using model (1) by an Interior Point method, knowing which k-beams will be used. This is the means to evaluate fitness, because better individuals are those which are capable of generating good dose distribution. The fitness is given by the inverse of the objective function value. A penalty is applied to the infeasible individuals, which hence receive a small fitness. (ii) Selection: The selection procedure used is roulette-wheel selection ( [13], [14]). (iii) Recombination: One position is randomly chosen and it represents the partition in the two parent individuals that will be divided to create the two new offspring. One of the offspring has the first part of the first parent, and second part of the second parent and the other offspring has the first part of the second parent and the second part of the first parent, as showed in the Figure 4. It is possible that the created offspring is infeasible, this occurs when all the genes are zero or it has more than desired number of ones. In this case, the infeasible chromosome receives a low value of fitness. If in this step the created chromosome is equal to another, it suffers a migration process, which means that it was thrown away and another feasible chromosome will take its place.  Then, the population is replaced and the steps are repeated until a specified number of generations is reached. The experimental results using this methodology will be presented in section 5.2.

Experimental results -a case study
The case study used is the cancer of the prostate region, as shown in Figure 2, using one slice of CT. After pre-processing, the image has 11,319 pixels. This image is obtained in Erasto Gaertner Hospital (Curitiba/PR/Brazil), with the project number being 2042.
The experiments are undertaken for two different instances. The first has a set of 8 possible beam directions, with 10 sub-beams each, and up to 6 beam directions could be chosen. The second has a set of 8 possible beam directions, with 10 sub-beams each, and up to 4 beam directions could be chosen. The possible beams directions for the two instances are shown in Figure 5.
In this configuration, the model (1) has 22,718 variables, and 11,319 constraints. Two different solution approaches are used. The first one is the exact methodology, purposed in [12], that will be presented in section 5.1. The exact solution results will be compared to the results obtained by the proposed matheuristic, which will be presented in section 5.2.
The software system Matlab was used, and run on an Intel PC with Core i5 CPU at 1.7 GHz.

Results of exact methodology
The multiobjective model (1)-(3) is solved with Branch and Bound, using the weighted sum method for the goals, represented on equation (4), with 0.4 as weight to the objective functions (c) and (d), ie γ 1 = γ 2 = 0.4, and 0.1 as weight to the objective functions (a) and (b), ie α = β = 0.1. These weights are chosen based on the experiments done in [12]. Two instances are proposed, using 8 possible beams to select up to 6, for the first instance, and 8 possible beams to select up to 4, for the second instance. The solution of these instances are shown on Table 1.  As presented by Table 1, this methodology is expensive in computational time. The matheuristic proposed which is a hybridisation of GA and Interior Point is used to compare the optimum result obtained with the exact method, and the experimental results are shown in the next section.

Results of matheuristic
The matheuristic described in section 4 is used for the same instances proposed in section 5.1. For each instance described in Table 2, we run 5 times the experiments, in which the process of selection selects 80% of the individuals, and the mutation probability is 5%.  In the fitness evaluation, the model (1) is used in the same scenario as in 5.1, with weights: γ 1 = γ 2 = 0.4 and α = β = 0.1. To solve it, the Interior Point methodology is used, and the results obtained are shown on Table 3. Table 3 presents the average results for each instance over the five runs, in terms of: objective function for each deviation dose, the weighted sum of the functions (represented by g), the average computation time, the percentage improvement in computation time relative to the exact method and the average number of generations to convergence.  Table 3. Average results of GA   The results on Table 3 show that the matheuristic reaches the optimum value in 20% of the experiments, for the first instance, and 8%, for the second one. Table 4 presents the best solution of each experiment in each instance. In the first instance, just one experiment do not reach the optimum, but has a quite small deviation to it (0.01%). In the second instance, for the five experiments done, three do not achieve the optimum, and have, on average, 0.9% deviation.
Both instances present better results using 5 generations and 15 population size, with 0.28% average deviation to the optimum, in the first instance, and 0.85%, in the second. The worst result has average 5.26% deviation to the optimum, in the first instance, and 3.22%, in the second.
The improvement in terms of computational time is very substantial, considering that it is between 80.34% and 95.49%, in the first instance, and between 71.88% and 94.98%, in the second, as presented on Table 3.
Even with a small number of generation, the results are impressive, whereas the experiments that converge have done it in less than 6 generations.
The population size is important to the convergence process, because better results are found with its increase, but have the disadvantage of also increase the computational time.
The presented instances shows the quality of the matheuristic to solve the combined beams selection and dose distribution problem in radiotherapy planning. It is faster than the exact method and finds the global optimum or very close to it in many executions done. More experiments in terms of crossover and mutation process could be conducted in order to improve the GA in this case study. Accordingly, other instances could be tried.

Conclusion
The planning of radiotherapy is a process that could be improved with operational research techniques. This paper introduces a methodology to find the beam directions for a set of possible beams to be chosen, and the optimum distribution dose for the selected set of beams.
We have conducted experiments to evaluate the performance of the proposed matheuristic in terms of the solution quality and computation time, comparing it to the exact method. The matheuristic results obtained a great improvement in computation time, and good quality solutions.
Although the matheuristic has achieved a substantial improvement in computation time, we still believe that it could be further improved by using a single population metaheuristics. So, a future work is to investigate other single solution metaheuristics to solve the problems of beam directions and dose distribution, with the intent to improve the computation time of the matheuristic proposed in this paper.
Inclusion of non-coplanar beams with CT 3D configuration is also crucial to the possibility of the use of these methodologies in radiotherapy planning in hospitals.