A multiphase direct aperture optimization for inverse planning in radiotherapy

NUCLEAR PHYSICS AND INTERDISCIPLINARY RESEARCH

A multiphase direct aperture optimization for inverse planning in radiotherapy

WANG Jie，

PEI Xi，

CAO Rui-Fen，

HU Li-Qin，

WU Yi-Can

Nuclear Science and Techniques

Vol.26, No.1

Article number 010502

Published in print 20 Feb 2015

Available online 20 Feb 2015

DOI：10.13538/j.1001-8042/nst.26.010502

57801

Optimization of the inverse planning becomes critical because it follows the invention of intensity modulated radiotherapy (IMRT) to shorten the previous "trial-and-error" treatment process and increase efficiency. In this paper, the inverse planning is used to direct aperture optimization in the ARTS (Accurate/Advanced Radiotherapy System). The objective function was quadratic, both tolerance and dose-volume constraint types are supported. The memory efficient conjugate gradient algorithm is used to cope with its large data. Furthermore, to fully exploit the solution space, a shortest path sub-procedure is coupled into the whole algorithm, thus giving further possibility decreasing the objective function. Two clinical cases are tested, indicating that the applicability of this algorithm is promising to clinical usage.

RadiotherapyInverse planningConjugate gradientShortest path

I. INTRODUCTION

Radiotherapy as an effective tool to cure cancer has been practiced for decades. Its guiding principle is to give sufficient uniform dose to a tumor target while sparing the normal organs. In the old days, this was not an easy task, and the treatment was basically based on the physician’s experience. The turning point occurred at the 1980’s when Brahme solved a simplified inverse problem analytically [1], and an ideal goal could be achieved with an inverse treatment to the radiotherapy problem. After that, different models and methods were used to investigate this issue. Another milestone is the introduction of "direct aperture optimization" [2], which simplified the traditional optimization process yet still guarantee a promising result. This work falls into the category of direct aperture optimization. The choices of objective functions are not limited. Quadratic deviation is a natural choice and was adopted by many researchers. Other popular objectives include linear [3], biological-based [4], and mix integer. A natural choice to optimize would be heuristic methods. Simulated annealing was adopted first in the optimization of radiotherapy by Webb [5], and other methods like genetic algorithm, evolutionary algorithm etc. are abundant. The natural advantages of heuristic methods are their ability to handle complex, non-regular objective functions, as they do not need gradient information to perform optimization. The shortcoming is that the optimization takes much time. While deterministic methods surpass in this aspect, by using gradient information of objective functions and careful implementation, they usually have a much shorter optimization time. Theories on deterministic optimization are enormous from which effective methods can be adopted on the problem of radiotherapy optimization. Due to the fact that many optimization algorithms only work on small-scale or middle-scale problems, they do not perform well on large-scale problems, and the computer would not even store all the data needed in the algorithms. The most effective method is conjugate gradient (CG) method to this kind of problem. Only part of the data needs to be stored in optimization, and it converges in a determined number of iterations. But still, special modifications shall be made to customize the standard conjugate gradient algorithm applicable to this problem. A sub-procedure based on shortest path problem is embedded into the optimization to improve the outcome. It searches the new apertures which have the potential to decrease the objective function value further more.

II. MATERIAL AND METHODS

A. Problem formulation, notation

In the Fusion Driven sub-critical System team [6-9], a radiotherapy system called ARTS [10-22](Advanced/Accurate Radiotherapy System) is developed. The dose engine, MCFSPB (Monte Carlo finite size pencil beam) [19], generates the dose kernel matrix for inverse planning optimization. We consider two types of organs, tumor target (T) and normal organ (NO). Sample dose are picked to fully represent the whole organ. Beam intensity, which is denoted by X, is a vector with thousands of components, and dose kernel is denoted by D. Beam direction is {B}. Each beam direction is decomposed into several small beam fields, which are call apertures, and denoted by {A}. Each aperture corresponds to an aperture weight, ω, and is described by the multi-leaf collimator’s leaf position, Y vector. A quadratic based penalty function was adopted, as follows:

f (d) = δ_{T} ρ_{T} (d_{T} - P_{T}) (d_{T} - P_{T})^{2} + δ_{NO} ρ_{NO} (d_{NO} - P_{NO}) (d_{NO} - P_{NO})^{2},

(1)

d = D x,

(2)

x = \sum_{i = 1}^{B} w_{i} \sum_{j = 1}^{A_{i}} ϕ_{j} (y),

(3)

δ (x) = {\begin{matrix} 1 x \geq 0 \\ 0 x < 0 \end{matrix},

(4)

ω, x \geq 0.

(5)

δ(x), which is a Heaviside function as an indicator, shows wether the sample point should be given to penalty. Meeting the dose prescription at the sample point is not penalized, and δ(x) =0, or else, penalized, and δ(x) =1. ρ is the radio-sensitivity factor of different organs. P represents physician’s prescription dose to the organ, and means upper tolerance, low tolerance, or mean dose in different context.

In Eq. (3), all the apertures gives the beam intensity vector involved in Eq. (1). The dimension of dose kernel is typically around 5 000×10 000, and poses limit to even the computational power of a modern PC.

B. Constraints

Optimization problems and radiotherapy inverse planning involves constraints. A group of normal organs called OAR (Organ at Risk) are supposed to be spared as much as possible, and various constraints are imposed, and an upper tolerance dose is limited by

D_{OAR} x \leq P_{OAR} .

(6)

If the geometrical configuration is too complex, the constraints would degrade the treatment outcome in tumor target, the physicians would like to weaken the constraints on OAR’s. This option is called as dose-volume constraint, including the dose d, and the volume percentage ρ%. For a typical dose-volume constraint, the percent volume of the organ receiving a dose over d should be under ρ%. In the objective function, the target volume received a dose d should be at least ρ%. Thus, the dose-volume constraint is more flexible than the tolerance limit. Small fraction of the OAR to be overdosed may lead to a uniformly high dose distribution.

The nonlinear dose-volume constraint incorporating into optimization will complicate the optimization process. Here, a straightforward sorting-based treatment of dose-volume was used, the underdose or overdose part will be penalized after sorting the dose.

C. Algorithms

To alleviate the computational burden, a multiphase strategy is applied. Each beam direction is initialized with one aperture, the aperture weights are optimized by the Fletch-Reeves conjugate gradient method. When the first phase stops, the second phase of new aperture generating starts. New apertures are generated by solving a shortest path problem, which will be explained in Section 2.5, added into the aperture pool, and the whole optimization process continues. After generating 5 new apertures or so, the satisfactory optimization process is obtained. Fig. 1 shows flow chart of the optimization process.

Fig. 1.

Flow chart of optimization.

D. Aperture weight optimization

Simultaneous optimization of aperture shape and aperture weight is possible, Webb [5] used simulated annealing. Genetic algorithm and evolutionary algorithm are adoptable, the optimization takes longer time.

Coupling the aperture weight and aperture shape is the main reason why this problem is hard.

\begin{array}{l} \min f (x); \\ x =_{\sum}^{B}_{\sum}^{A} ω ϕ (y) . \end{array}

(7)

Handling them separately shall be a wise option. Fixing y makes the problem a standard quadratic programming problem, which is well solved. After aperture weight optimization, a shortest path problem is formulated as a means to further explore the solution space, which is confined by the previous limitation on aperture weight.

Conjugate gradient method is suitable for solving large sparse problem, and saves both memory and time. Normally, conjugate gradient method performs very well for a standard positive definite matrix based quadratic programming.

\min x^{T} H x + P x .

(8)

This problem is not quite standard. The objective function f(d) = δ_T(d_T - PT)(d_T - P_T)² + δ_NO(d_NO - P_NO)(d_NO - P_NO)² is a weighted quadratic form. Although the square term is positive semi-definite, the value of Heaviside function ahead each term is not alike each time. So, the problem is evolving with time. Nevertheless, the direction provided by a conjugate gradient algorithm is still useful for the evolving. The non-negativity constraint on x is incorporated in the optimization by setting the negative components to zero. The following conjugate gradient algorithm is used in this work:

(1) Set initial x₀, y₀ and tolerance ϵ;

(2) Calculate $\nabla f (x_{0})$ , if $| \nabla f (x_{0}) | \leq ϵ$ , stop; otherwise, go to Step 3;

(3) Let $p_{0} = - \nabla f (x_{0})$ , let k=0;

(4) Line search for tk, x_k+1 = xk + tkpk, if xj <0, set xj =0;

(5) Calculate $\nabla f (x_{k + 1})$ , if $| \nabla f (x_{k + 1}) | \leq ϵ$ , stop; otherwise, go to Step 6;

(6) $p_{k + 1} = - \nabla f (x_{k + 1}) + λ_{k} p_{k}, λ_{k} = \frac{| \nabla f (x_{k + 1} {) |}^{2}}{| \nabla f (x_{k} {) |}^{2}}, k = k + 1$ , go to Step 4;

(7) Stop.

Both tolerance constraint and dose-volume constraint are supported and handled by incorporating them into the objective function using Heaviside function. Dose-volume constraint is realized by penalizing the region overdosed or underdosed after sorting.

\begin{array}{l} g (d) = f (d) + δ_{tol} (d_{tol} - P_{tol}) (d_{tol} - P_{tol})^{2} + δ_{dv} (d_{dv} - P_{dv}) (d_{dv} - P_{dv})^{2}; \\ d = D x . \end{array}

(9)

E. Shortest path problem

A necessary condition of the minimization of the objective function is

\partial f / \partial x = 0.

(10)

The results obtained from Section II D do not meet the condition because compromises are made. To further explore the solution space, finding the solution potential missed can be achieved by adding the appropriate aperture grids into the solution. The incorporation of a sub-procedure is based on an intuition that the negative partial gradient of the objective function to the beamlet grid may lead to a decrease of the objective function. Gradient of the objective function to the beamlet, ∂f/∂x_i, is calculated, and the shortest path problem is used for each row of each beam field.

\min sum (\partial f / \partial x), - x \in B .

(11)

The connectivity requirement of aperture is directly constructed to the optimization process. Various strategies can be used to define what does shortest path mean. Here, a naïve strategy chosen is the smallest connected beamlet array has the most potential to decrease the objective function.The solution of this problem consists of the new aperture which will be added to the aperture pool for optimization. This sub-procedure is an effective means to further explore the solution space.

III. RESULTS AND DISCUSSION

The test cases consist of a prostate cancer case and a lung cancer case. Both of them have voxel sizes of (5 mm×5 mm×5 mm). The dose kernel was computed based on an Elekta linac 6 MeV electron beam profile in ARTS. The configuration of the computer used is: Intel Core2 Quad CPU Q9550 @2.83 GHz, 3.50 GB memory.

A. Prostate cancer

The first clinical case is a prostate cancer. It consists of 35 CT scans. Organs delineated are tumor target, rectum, bladder, and femoral head. Fig. 2 shows the distribution. Tumor target is in the middle of four normal organs. Bladder and rectum are especially close to the tumor target.

Fig. 2.

(Color online) Organ distribution of the prostate cancer.

Constraints of the organs are given in Table 1.

Constraints of prostate cancer

Organ	Prescription (cGy)
Tumor	5100 ࣘ p ࣘ 5500, mean dose 5200
Rectum	At most 30%, larger than 2000
Bladder	At most 30%, larger than 2000
Femoral head left	At most 10%, larger than 3500
Femoral head right	At most 10%, larger than 3500

The 7 beam fields are used in the case with an equi-spaced angle of 50, the total maximum apertures are 50. The optimization took 67 s. The iso-dose contour and DVH (dose-volume histogram) graph are given in Fig. 3.

Fig. 3.

(Color online) ISO-dose graph (a) and DVH graph (b) of the prostate cancer.

Although the tumor target is very irregular, the 100% curve encompasses the region very well. The 90% and 80% curves expanded the central region a little, but they avoided the rectum as much as possible: almost all the rectum region was avoided by the high dose field, only a tiny bit was exposed to it. Other organs also receive appropriate dose levels. In the DVH graph, the curve of tumor target is steep, indicating good dose conformity to the tumor target. The left and right femoral head meet the prescription. And rectum and bladder were constrained very well. The optimization outcome is satisfactory.

B. Lung cancer

The clinical case of lung cancer contains 53 CT scans. The organs delineated are tumor target, lung, heart and spinal cord. The distribution is shown in Fig. 4. The tumor target is in the left lung.

Fig. 4.

Organ distribution of the lung cancer.

Objective and constraints of the organs are given in Table 2.

Constraints of lung cancer

Organ	Prescription (cGy)
Tumor target	7000 ࣘ p ࣘ 7500, mean dose 7200
Lung left	At most 30%, larger than 3000
Lung right	At most 30%, larger than 3000
Heart	p ࣘ 4500
Spinal cord	p ࣘ 4700

There are 9 beam fields in this case with an equi-spaced angle of 40, the total maximum apertures are 50. The optimization took 376 s due to its larger size. The iso-dose contour and DVH graph are shown in Fig. 5.

Fig. 5.

(Color online) ISO-dose graph (a) and DVH graph(b) of the lung cancer.

The 100% iso-dose curve followed the same track as the tumor target contour, though its shape is not quite regular. The 90% and 80% curves expanded a little, but avoid the neighbor spinal cord completely, which is exactly what intensity modulated radiotherapy should achieve. An interesting fact is that the 80% iso-dose curve is elongated at the bottom direction because the tumor target sits in the bottom region of the left lung, and higher beam energy was given in the reverse direction resulting in a higher tail dose deposition. Other organs also received appropriate dose level. From the DVH graph, a steep curve of tumor target is presented, indicating good conformity, spinal cord, heart and lung left were confined very well under the dose constraints. Dose of lung right is not as low as the other three, but it still meets the dose criteria.

IV. CONCLUSION

Ever since inverse planning concept was introduced into the radiotherapy field, optimization seemed to be a necessary issue which needs to be addressed. Finding an acceptable treatment plan in an acceptable time became the main purpose of the radiotherapy treatment planning optimization. Unfortunately, the same as other computation issues in computer science, the two aspects are essentially in conflict to each other, thus need tradeoffs. Generally, optimization problems involved in radiotherapy are broad in type, large-scale and often hard to solve. This work addressed the issue as a direct aperture optimization using a deterministic method. Conjugate gradient algorithm adopted in this work enables us to solve large-scale problem efficiently. The result is promising, and can be obtained in a relatively short time. The heuristic part of the main algorithm may be the hidden reason why this algorithm performs better than some results of the previous literature. Some other forms of objectives and constraints can be explored in future, which could potentially have a better performance.

References

[1]

A Brahme.

Optimization of stationary and moving beam radiation therapy techniques

. Radiother Oncol, 1988, 12: 129-140. DOI: 10.1016/0167-8140(88)90167-3