Mixed Integer Programming Approach In Radiation Treatment Planning

Mixed Integer Programming Approach In Radiation Treatment Planning

S. P. Kishore 1, N. V. S. L. Narasimham 2

1(H & S Department,Vardhaman College of Engineering,Shamshabad,Hyderabad-501218,India)

2(Department of Mathematics,GNITS,Shaikpet,Hyderabad-500008,India)

ABSTRACT

Radiation is one of the most common treatments for cancer. It is often a part of the main treatment for specific types of cancer, such as, lung, head and neck, Hodgkin disease. Radiation can sometimes be used with the treatments, such as chemotherapy or surgery. We present mixed integer programming model of inverse planning for IMRT. The challenge in therapy planning is to determine a treatment plan which would involve finding irradiation directions, beam parameters and optimize intensity profiles. Radiation uses specific equipment to deliver lethal doses of radiation to the cancer cells. In the process of radiation, normal cells do get affected, but most of them fully recover in due course of time and go back to the working way. IMRT uses computer-generated images to plan and then deliver even more tightly focused radiation beams to cancerous tumors than is possible with conventional radiotherapy. Inverse planning in radiation therapy is a trial and error method with an ability to improve the outcomes of treatment to cancer patients. The aim of optimization technique is to deliver a dose of radiation to the cancer cells as possible, with minimal damage to the surrounding vulnerable normal tissues. The proposed model has an advantage over the existing models and an attempt has been made to improve the quality of the treatment plan within a realistic time frame.

Keywords-intensity modulated radiotherapy (IMRT); inverse planning; Multileaf collimator (MLC); mixed integer linear program (MILP);beam angle selection;genetic algorithm(GA).

1. Introduction

   Traditional radiotherapy uses forward approach and inverse planning. In the forward approach, beam directions and intensity profiles are selected by the radiation therapy planner and the dose distribution is calculated. In addition, the chosen intensity profiles should be satisfactory or else the process is repeated until a desired dose distribution is obtained. Inverse planning is a backward approach to radiation treatment planning which takes into consideration the target dose as the starting point and search for the intensity profile that will produce that dose. Inverse planning does not make any prior assumptions which make it the best method to produce more improved treatment plans. Inverse planning clearly specifies the objective of the problem, it specifies a desirable dose distribution in the target and an acceptable dose distribution in the critical structures

   IMRT treatment planning is performed in three phases. The first is the beam orientation selection which determines the selection of the gantry positions and their angles. The second part of the process is called fluence map optimizations which discretise each beam into small bixels and determine the intensity profile for each beam that yields the best quality treatment plan. The third part of the process is MLC leaf sequencing which is needed to

   decompose the intensity profiles of each beam into a collection of deliverable apertures and corresponding intensities. The application of a MultiLeaf Collimator (MLC) has enabled the use of complex intensity profiles for a field. The MLC consists of narrow blocking leaves which can be moved forwards and backwards in a treatment field (see[2])

2. Mathematical model

   In order to present the mathematical model, we need to present some fundamental assumptions. We need to discritise the three dimensional volume of the patients body containing the tumor, nearby organs at risk and any other parts of the body where the beams might enter or exit through. This volume is descritised into equidistant volume elements voxels and the beam elements bixels.(see Hamacher et al.[1]) In clinical practice MRI or CT scans of the part of the body are used to create 3D simulation of the body.Let us assume that K body parts are under consideration where k=1 indicate the tumor or the targer structure and the indices k=2, 3, 4..K indicate organs at risk.Typically the order of K is between 3 and 7.Let M=M1 +M2 +……..+Mk be the number of voxels where each voxel is uniquely identified as belonging to target structure or normal tissue .During the treatment, each voxel will absorb a dose of radiation. We denote the dose distribution in the body by a K dimensional vector DR K .

   The Standard IMRT model is given by D = AX

   or

   n

   n

   Di aij x j ,i=1,2,..k

   j1

   Where A is the matrix such that

   aij

   represents the amount of dose in voxel i per unit intensity from the jth bixel.

   ARk x n is known as the influence matrix.

   This linear dose model D=AX is an approximation of the true nature of the radiation treatment.The matrix A is computed by using different approximation methods.Dose calculation is done using many faster,but less accurate methods and Monte carlo sampling techniques(see Francescon P et al. [3]) are among the popular one.we assume that the matrix A has non zero elements,so all the voxels will receive some radiation and every bixel influences one of the voxels dose.

   1. Mixed integer programming

      The Most favorable approach to IMRT is to express gantry positions in the optimization model and to solve it .This is supposed to be a relatively simple approach to formulate and mention all the aspects of the model being developed. In this paper,we have used binary variables to control the number of irradiation directions in the solution,which resulted in a mixed integer linear program(MILP).In MILP formulation we have used both continuous and Binary variables.In binary problems,each variable can take on the value of 0 or 1.This represents the acceptance or rejection of an option,the turning on or off the various gantry angles.Radiologist generally use a cumulative dose volume histogram(DVH) to determine the quality of a radiation treatment plan.Figure 2 illustrates the cumulative DVH in which the fraction of the patient that receives at least a specified dose level

      K

      K

      The mathematical model is as follows Min F(x,T)=wk Tk ,

      k=1

      where the parameters wk , k=1K indicate relative importance of tumor and organs at risk

      In this formulation the constaints (1) represent calculation of dose distribution in the tumour and equires that the dose is greater than or equal to the desired dose with some devaition T1 ,where I is a vector of all ones with an appropriate dimension.Similarly (2) requires dose in organs at risk(OAR) to be less than or equal to upperbound with deviation Tk,L1 represent lower bound on the tumour and Uk is an upper bound for OAR.T=(T1,T2,..,Tk) measures the deviation between delivered dose and the desired dose levels in each volume of interest.The constraints

      (3) ensures that the values are non-negative.

      D1 =A1X (L1-T1) I , (1)

      Dk =Ak X (Uk +Tk ) I

      (2)

      All

      x 0, Tk 0

      (3)

      A mixed integer formulation (see [5,6,7,11]) is brought by introducing the Boolean variables yh which indicate the acceptance and rejection of an option indicated by (4). The constraints (5) and (6) are introduced to control the number of irradiation directions and beam angle selection is done accordingly so that the radiation is emitted fom selected directions.

      y = 1;radiation is emitted

      h 0;otherwise

      (4)

      <>H

      H

      yh R

      h=1

      xhi yh ,i 1, 2,…..N; h 1, 2…..H

      (5)

      (6)

      Let us now investigate the size of the MILP formulation.Let us consider a simple example by taking 72 possible beam orinetations descritised with 5o intervals which would result in 72 boolean variables yh denoting which of the beams would be turned on.For each beam varaible we have around 150 beamlet variables,xi making it a total of 11,800 beamlet variables.In addition to that the formulation take into account around 100000 voxel variables and 40000 constraints also.This was put on the state of art commercial solver CPLEX 8.0 for solution and it took almost 3 days to arrive at a solution.It appears that MILP of this type are difficult to solve by general purpose codes so we limit beams angle selection toa minimum.

   2. Dose-Volume Constraints

      Mostly all the recent IMRT inverse planning systems allow the specification of dose-volume constraints (DVCs). Typical dose-volume constraints limit the relative volume of a structure that receives more or less than a particular threshold. The cumulative dose volume histogram displays the percentage of volume of the patient that receives a relative dose (see figure 2).In few cases the radiation oncologist does sacrifice a portion of a region at risk in order to improve the chances of curing the disease. For example, the constraints can be specified of the form No more than x % of OAR can exceed a relative dose of y.This kind of requirement is very useful for the clinic and is termed as partial volume constraint.

      Figure 2.Cumulative dose volume histogram

   3. Genetic Algorithm

      Genetic alorithm (GA) work well on mixed (continuous and discrete) combinatorial problems.They do not stuck at local optima unlike the gradiant search methods((see, e.g., Erhgott et al. [2]).GA are the heuristic search and optimization techniques based on the principles of natural selection and evolution.In GA algorithm mutation technique is used to obtain new generation of samples. Mutation is done by Gaussian probability distribution with

      zero mean N(0,i ) The GA algorithm is illustrated below

3. Discussions

   All commercially available IMRT treatment planning systems, as well as most of the research in the medical physics literature to date, use local search (such as the conjugate gradient method; see, e.g., Shepard et al. [6], Xing and to find a satisfactory treatment plan. In a recent comprehensive review of the radiation therapy literature, Shepard et al.

   [6] surveyed many techniques previously used, including some simple linear and convex programming formulations of the problem. More recent operations research approaches are multi-criteria optimization (see,e.g., Hamacher and Kufer [1]; Bortfeld et al. [8], and mixed-integer linear programming (see, e.g., Erhgott et al. [2]).GA is the most effective method till date of searching the sets of gantry angles.It was simultaneaously run on k machines to arrive at the next generation of solutions,so as to reduce the time constraint

4. Conclusions

In clinic use , the common practice is to use more beams in the treatment ,but that results in higher verification and treatment times to avoid that, we keep the number of beams used to be as low as possible. In this mathematical model the focus is on the important aspects of inverse planning unlike the trial-and error approach. Optimization of the beam angles resulted in an improved treatment for IMRT.The MILP was solved using the state of art commercial solver CPLEX 8.0 making it as a powerful tool for radiotherapy planners. We sole the LP relaxation of the mixed integer program to obtain an optimal solution. Based on this model, we had overcome some of the problems faced with optimization techniques in IMRT treatment planning like the use of priority weights, beam intensity profiles, overdosing the organs at risk (OAR), under dosing the cancerous cells. GA is one of the popular methods to find the space of gantry angles .A number of challenges still need to be investigated. Model formulation and implementation of different optimization techniques is still required to solve problems of this size within a realistic time frame.

References

1. Hamacher, H.W. and Küfer, K.H. (2002): Inverse Radiation Therapy Planning a multiple objective optimisation approach, Discrete Applied Mathematics, 118:145-161.

2. Ehrgott, M. and Burjony, M. (2001): Radiation Therapy Planning by Multicriteria Optimisation, Proceedings of the 36th Annual ORSNZ Conference, 244-253.

3. Francescon P, C ora S, C hiovati P. Dose verification of an IMRT treatment planning system with the BEAM EGS4-based Monte C arlo code. Med Phys. 2003;30(2):14457

4. Zhang X, Liu H, Wang X, Dong L, Wu Q, Mohan R. Speed and convergence properties of gradient algorithms for optimization of IMRT. Med Phys. 2004;31(5):114152.

5. E. Lee, T. Fox, I. Crocker, Integer programming applied to intensity-modulated radiation treatment planning optimization,Ann. Oper. Res. 119 (2003) 165181.

6. D. Shepard, M. Ferris, G. Olivera, T. Mackie, Optimizing the delivery of radiation therapy to cancer patients, SIAM Rev. 41 (1999) 721744.

7. Hager, S. Huang, P. Pardalos, O. Prokopyev (Eds.), Multiscale Optimization Methods and Applications, Springer,2005, pp. 205228.

8. T. Bortfeld and W. Schlegel, Optimization of beam orientations in radiation therapy: Some theoretical considerations, Phys. Med. Biol., 38 (1993), pp. 291304.

9. Boland N, Hamacher H W, and Lenzen F. 2004. Minimizing beam-on time in cancer radiation treatment using multileaf collimators. Networks 43(4):226- 240

10. Y. Zhang, M. Merritt, Dose volume-based IMRT fluence optimization:A fast least squares approach with differentiability, elsevier,2008, pp. 13651387.

11. Lenrick A.Johnson,Optimization of Irradiation directions in IMRT Treatment,University of Auckland,2004.

12. Steuer, R. Multiple Criteria Optimization: Theory, Computation, and Application, Wiley 1986