Development of a test bench for Gamma Knife Optimization

(1)

Development of a test bench for

Gamma Knife Optimization

V I C T O R I A S V E D B E R G

(2)

(3)

Development of a test bench for

Gamma Knife Optimization

V I C T O R I A S V E D B E R G

Master’s Thesis in Optimization and Systems Theory (30 ECTS credits) Master Programme in Mathematics (120 credits) Royal Institute of Technology year 2014

Supervisor at KTH was Johan Karlsson Examiner was Johan Karlsson

TRITA-MAT-E 2014:02 ISRN-KTH/MAT/E--14/02--SE

Royal Institute of Technology

School of Engineering Sciences

KTH SCI

(4)

(5)

Abstract

The dose distribution in the Gamma Knife is optimized over the weights (or Beam-on time) using different models other than the radiosurgical one used in Leksell Gamma Plan R

(6)

(7)

CONTENTS CONTENTS

1 Introduction

The Leksell Gamma Knife (developed and produced by Elekta) is an instrument for treating braintu-mors and various malfunctions in the brain. To optimize the given dose in a treatment using the Leksell Gamma Knife, a radiosurgical model is used. In radiosurgery, high-energy radiation is aimed at the target (most often a tumour) in order to destroy the target cells. The radiosurgical model is based on coverage, selectivity and gradient index to achieve the aim, which is to cover the target (and only the target) with the planned dose and make the dose drop rapidly outside the target for it to quickly drop, thus harming the healthy tissue outside the target as little as possible. In this model no account is taken for the riskstructures in the brain and if the target is big, the planning-time tends to become large due to the large number of degrees of freedom.

The goal when using radiotherapy when treating a tumour is to destroy the DNA, which In radiotherapy, the models which are based on the Dose Volume Histograms (DVH), are the most common for optimizing a treatment plan. DVH describes the percentage of a target or a risk-structure which achieves a certain dose. Because of the non-convexity of the DVH and in order to use fast optimization techniques, sim-plifications have to be done. These simsim-plifications are called DVH-surrogates and are widely used in the field of radiotherapy, e.g. at Princess Margaret Hospital (PMH) in Toronto, researchers have replaced the radiosurgical treatment plans in the Gamma Knife with DVH-surrogates.

Radiobiology aims to explain the biological response in the tumor tissue and in the organs at risk when irradiatated with high-intensity ionizing radiation [12]. Radiobiological models are more and more com-mon in cancertreatment devices. These models are usually based on Tumour Control Probability (TCP), Normal Tissue Complication Probability (NTCP) and Equivalent Uniform Dose (EUD). One problem using radiobiological models is the large uncertainties in the parameters. Along with the number of parameters used these models may become very uncertain in predicting the outcome of the treatment. In some tumours there may arise areas where there are less molecular oxygen than in the rest of the tissue. These areas are called hypoxic areas and are more resistant to radiation and higher dose must be given in these areas to achieve tumour control. The model parameters for TCP change for it to express the stronger radio resistance.

(10)

2 BACKGROUND

2 Background

In this thesis, models based on Radiobiology and Radiotherapy are tested in the Radiosurgical instrument, the Gamma Knife. Therefore, a brief introduction into Radiosurgery, Radiotherapy and Radiobiology will be given and along witg a technical overview of the Gamma Knife.

2.1 Radiosurgery

B¨orje Larsson (professor in Radiobiology at Uppsala University) and Lars Leksell (neurosurgeon at Karolinska University Hospital) developed the concept of radiosurgery in 1951 [19]. This concept was to direct a single high dose of radiation to the intracranial region of interest stereotacticly. The word stereotactic is commonly used as a prefix to radiosurgery and refers to the three-dimensional coordinate system which makes it possible for the treatment planner to identify and use the coordinate systems of the diagnostic images and the actual position of the patient, thus simplifying the treatment and raising the accuracy of hitting the planned shot position (also called isocenter). The treatment plan is done in the Leksell Gamma Plan.

The fundamental principle of radiosurgery is that of ionization of target tissue (for instance a tumour), by means of high-energy radiation, usually gamma radiation. By ionization the tissue, the number of ions and free radicals (an atom, molecule or ion with unpaired valence electrons) which are harmful to the cells increase (more about this in Section 2.3.1). These ions and radicals can produce irreparable damage to DNA, proteins, and lipids, resulting in the cell’s death. Thus, biological inactivation is carried out in a volume of the treated tissue. The radiation dose is usually measured in Grays, where one Gray (Gy) is the absorption of one joule per kilogram of mass.

Radiosurgery is performed by a team of neurosurgeons, radiation oncologists, and medical physicists to operate and maintain very sophisticated, accurate and complex instruments. The precise irradiation of targets within the brain (and upper part of the spine) is planned using information from diagnistic images which are obtained via computed tomography, magnetic resonance imaging, positron emisson tomography and angiography. An advantage of stereotactic treatments is the delivery of the right amount of radiation to the tumour in one or a few fractions compared to traditional treatments, which is often delivered in 30 fractions over a timeperiod up to 10 weeks. Also, treatments are given with extreme high accuracy, which limits the effect of the radiation on the surrounding healthy tissues due to the sharp dose gradients. One problem with stereotactic treatments is that they are only suitable for small and medium sized tumors, since the treatment time may increase considerably for large tumours.

2.1.1 Leksell Gamma Knife

The Leksell Gamma Knife (or shortly the Gamma Knife) is a medical device for treatment of brain tumours, developed by the company Elekta AB. The Gamma Knife was first developed by B¨orje Larsson, Ladislau Steiner and Lars Leksell in 1967 and is based on the concept of radiosurgery. The most recent model is Leksell Gamma Knife PERFEXION R _{presented in Figure 1. It was introduced in 2006. in this} section a small technical overview of the Gamma Knife is given.

(11)

2.1 Radiosurgery 2 BACKGROUND The most recent model is Leksell Gamma Knife PERFEXION R _{presented in Figure 1. It was introduced} in 2006. The Gamma Knife treats benign and malignant tumours and malfunctions in the brain by aiming gamma radiation on the tumour cells. The radiation source in the Gamma Knife PERFEXION R is Cobalt-60 and each of the sources is distributed in five circular rings in the device which is divided into 8 sectors. By radiating from each sector the gamma radiation is aimed at a target in the brain, hence focusing the rays to a so called isocenter as in Figure 2. The isocenters may have different shapes. The Gamma Knife is heavily shielded by lead, iron and tungsten to avoid leakage of gamma radiation. To avoid errors from patient movement during the treatment an aluminum frame is surgically fixed to the skull. This frame is fixated to the Gamma Knife in order for the head to remain stationary during the treatment. The frame becomes a reference for the coordinate system when calculating dose (hence a Stereotactic frame).

Figure 2: Conceptual picture of the Gamma Knife.

Focusing the radiation to an isocenter is the most important objective of the Gamma Knife (and in radiosurgery overall). By focusing the radiation, the intensity outside the focus point can be low enough to spare the tissue from any severe damage, but in the focus point become high enough to kill the tumour cells. Hence, the tumour cells die while the healthy tissue are relatively spared. This can result in shrinkage or a complete disappearance of the tumour. Gamma Knife radiosurgery has proven effective for patients with benign or malignant brain tumors up to 4 centimeters in size, vascular malformations such as an arteriovenous malformation (AVM), pain or other functional problems [11], [17], [14], [25]. The risks of gamma knife treatment are very low, and complications are most often related to the condition being treated.

2.1.2 The Gamma Knife Treatment Plan

(12)

2.1 Radiosurgery 2 BACKGROUND • Isocenter-packing

• Full optimization

Solving the isocenter-packing problem1 _{gives the positions of the isocenters (the coordinates where the}

gamma rays are converging as well as their shape) and hence the initial number of isocenters are given. The algorithm is constructed as follows: first the target is covered by as large shot as possible and the isocenters are placed so that it touches the target volume periphery, without overlapping other shots too much. When no more shot positions exists, even for the smallest volume, the volume covered so far is treated as non-target and the algorithm is repeated with the reduced target volume starting with the largest shot. Thus the target is filled from the surface and inwards [9]. A short graphical description of this algorithm i given in Figure 3. Only a subset of all isocenter shapes (isocenter states) are used and a shape of an isocenter is defined by {x|x ∈ R3, D(p, x) ≥ 0.5Dmax}, where p is the isocenter position, D

is the dose and depends on the size of the shot and Dmax is the maximal dose.

Figure 3: Illustration of the isocenter packing. A subset of all isocenter shots are used. The aim is to fill out the boundary as much as possible, then remove the isocenters, and fill the new volume. In the first (left) isocenter-packing, all isocenter shapes are allowed on the first boundary. In the right one the largest shot is not allowed on the first boundary. Image from [9]

Now the optimization problem may be solved. It is based on coverage, selecticity, gradient index and Beam-on time. These are stated in Equations 1, 2, 3 and 4. Here PIV is Planning Isodose Volume, the tissue covered by the planned isodose, and TV is Target Volume, i.e. the volume of the target (for example a tumour). ISO is the level of the isodose in percentage. Now, isocenter positions, isocenter-states and weights are varied freely.

Figure 4: Wenn-diagram of PTV and TV. Image from [9]

Coverage: C = V (P IV ∩ T V )

V (T V ) (1)

1_{More info}

(13)

2.2 Radiotherapy 2 BACKGROUND

Selectivity: S = V (P IV ∩ T V )

V (P IV ) (2)

Gradient index: GI = V (P IVISO/2) V (P IVISO)

(3)

Beam-on Time: Tbeam−on= Niso X

i=1

Tbeam−on,i (4)

CostF unc = C

min(2α,1)_{· S}min(2(1−α),1)_{+ βGrad + γT ime}

1 + β + γ

where α, β, γ ∈ [0, 1] are the weights befined by the user.

(5)

Here V (A) is the volume of the set A, Niso is the number of isocenters and Tbeam−on,i is the beam-on

time from isocenter i. The objective is demonstrated in Equation 5 where Grad is ae function of the gradient index and T ime is a function of the beam on times. Organ at risk are not being considered in the objective function. Instead poor selectivity is penalized so all tissue outside the target volume is treated as risk structures. A poor gradient index can also be penalized to create a steeper fall off of lower isodoses.

For the dose calculations two different calculations are performed. Those are • Convolution algorithm

• TMR-10 model

TMR-10 algorithm2 _{models all tissue in head as water [10]. This is proven to be good for targets in}

center of brain, but less accurate for heterogenities in the brain and peripheral regions. The convolution algorithm3_{takes tissue heterogenities into account in the dose calculation [8] and is only used to check the}

dose plan. However both give a mutually consistent result in a water phantom. During the optimization in the Gamma Knife, only the TMR-10 algorithm is used.

2.1.3 Leksell Gamma Plan R

The Leksell Gamma Plan R _{is the software for the treatment planning system for the Gamma Knife. It} provides a user friendly work system to the planner where all diagnostic images of the patient are collected and the dose delivery is planned. Using the CT-scan, the MR-image is fitted into the coordinates of the Gamma Knife and slices of the head are displayed. In these it’s possible to outline the tumour and organs at risk and perform an placement (organs at risk are only consider in the isocenter-placement). Then an optimization of the dose-distribution is performed, during which one may place out extra isocenters or remove an isocenter until a good plan is created. To investigate a plan isodoses and DVH:s are displayed. This plans are usually done by a neurosurgeon, or a medical physicist, but the plan must be approved by a neurosurgeon or an oncologist. When a plan is approved, the plan is sent to the Gamma Knife and executed.

2.2 Radiotherapy

Radiotherapy, similar to radiosurgery, also utilizes ionizing radiation to kill cells by focusing one or more beams to a planning target volume. In Conventional Radiationtherapy (CRT) the beams have a uniform intensity field which makes it harder to shape the field to avoid organs at risk (OARs). A progress was made in the 3D treatment planning by the introduction of Multileaf Collimators (MLC:s) to radiotherapy (first proposed by Brahme [3]). This was the beginning of Intensity-Modulated Radiotherapy (IMRT) where the radiation intensity across the beams can be modulated. This allows a better possibility of shaping the dose distribution of the planning target volume and sparing the OARs and normal tissue. The difference in dose shaping between CRT and IMRT is demonstrated in Figure 5.

2_{More info}

(14)

http://www.elekta.com/dms/elekta/elekta-assets/Elekta-Neuroscience/Gamma-Knife-Surgery/pdfs/LGP-2.2 Radiotherapy 2 BACKGROUND

Figure 5: Difference in dose shaping between CRT and IMRT. Image from IAEA.org

Figure 6: Illustration of MLC:s in operation.

Before radiation from a beam reaches the patient it has to pass through a MLC. The MLC consist of a number of pairs of metal leaves, usually made of tungsten, which can move and position itself into different shapes in order to block the radiation passing through the metal leaves and hence shape it. In Figure 6, the idea of dose shaping using multileaf colimators is presented. In order to model the beam intensities mathematically one discretize the beam into small ”bixels” or ”beamlets”.

There is a close resemblance between radiosurgery and radiotherapy but is in a fundamental plan not based on the same concepts. The aim of stereotactic radiosurgery is to kill target tissue with high dose while preserving the surrounding normal tissue utilizinf the sharp gradients of the Gamma Knife, where fractionated radiotherapy with less sharp gradients relies on a different sensitivity of radiation of the target and the surrounding normal tissue to the total accumulated radiation dose.

2.2.1 Dose-Volume Histogram

Dose Volume Histograms, abbreviated a DVH is a commonly used in clinical practice in Radiotherapy and Radiosurgery due to its economical and easy way to represent the entire dose distribution in a structure. The lost information is the spatial resolution of the dose distribution. It was first suggested by Bortfeld [1] and it is based on the physical dose. The planning structure considered (target or organ at risk) is divided into a number of voxels. The set of all voxels v in the planning structure S which attains dose d0 can be calculated using Equation 6.

vol : d07→ p(d0_{) :=} V ({v : dv = d0})

V (S) (6)

In other words the DVH provides the volume percentage of the planning structure in which at least the dose d is attained, i.e. the DV H(d) can calculated using Equation 7.

(15)

2.2 Radiotherapy 2 BACKGROUND The DVH is highly non-convex and non-convave. A simple example is that if x,y is the value of the dose in two voxels repsresenting the planning structure S. If x ≤ d and y ≤ d then 1₂(x + y) ≤ d. If x ≥ d and y ≥ d then 1₂(x + y) ≥ d. But what is the case of 1₂(x + y) ≥ d when x ≤ d and y ≥ d or the opposite x ≥ d and y ≤ d. Convex optimization is not applicable in this case. Because of this, transferring the DVH into more computationally desirable properties is an important field of ongoing research. There are some ways to control the DVH. Firstly, one simple possibility is to add a maximum or minimum dose constraint in the target and/or in the organs at risk. In Figure 7 a maximum and minimum constraint on the dose is demonstrated in a DVH-curve.

Figure 7: A maximum and minimum constraint on the DVH-curve. Image from [1]

For some organs maximum dose constraints is not meaningful at all. Organs like lungs and kidneys don’t have a large volume effect. If a part of these organs fail to function, the other part can still function without any major problem. In other words the tolerance dose in lungs and kidneys are not very dependent on the irradiated volume percentage. In the organs with a large volume effect, for instance the eye nerve, if a small piece of volume is irradiated, the tolerance dose is considerably lower than for larger irradiated volumes. In the organs with a large volume effect one can utilize the properties of Dose Volume Constraints, or DVC is also a common way to put a constraint on a DVH. The constraint is that no more than Vmax of the volume should receive more than a dose Dmax. In Figure 8 this is visualized

as a barrier with a corner at the point (Dmax, Vmax).

Figure 8: Visualization of a DV-constraint. Image from [1]

Unfortunately DVC-constraints lead to non-convex objective function [7][29]. Although it is shown that the resulting local minima in the overall objective function in simulations are of minor practical relevance [20].

(16)

2.3 Radiobiology 2 BACKGROUND

Figure 9: By using a small penalty factor on the maximum constraints some excess dose is allowed beyond the constraint. By a larger penalty factor prevents any larger dose than Dmax. Image from [1]

2.3 Radiobiology

The aim of radiobiology is to describe what happens when a cell is irradiated. Reasearch is made on this area today since it can contribute to the development of radiotherapy and radiosurgery in three important ways.

• It may provide an extended conceptual basis for radiotherapy and radiosurgery, by identifying mechanisms and processes that underlie the response for tumours and normal tissue to radiation and help to explain observed phenomena. Examples are knowledge about hypoxia, reoxygeniation and repair of DNA-damage (all these are explained later in Section 2.3.1).

• Development of treatment plans for radiotherapy and radiosurgery, for instance hypoxic cell sensi-tizers, to decrease the effect of an hypoxic environment (more in Section 2.3.5).

• Advice on the choice of schedules for radiotherapy and radiosurgery. For instance conversion formulas for changes in fractionation or dose rate. It may also provide a method for predicting the best treatment for the individual patient.

The role of oxygen is one positive example that has led to benefits in treatments today. The current method is empirical and further knowledge about cellular and molecular nature of radiation effects will lead to further development of the radiotherapy and radiosurgery [12].

2.3.1 What happens in an irradiated cell

When a cell is irradiated a number of processes commence. In Figure 10 the different phases and the time-period when they occur is graphically illustrated.

(17)

2.3 Radiobiology 2 BACKGROUND At first there is a physical phase where interactions occur between charged particles and the atoms of which the tissue is composed. A charged particle passes some living tissue and as it does, it interacts with the orbital electrons. Some is ejected from the atom (ionization) and some is raised to a higher energy level (excitation). If these are sufficiently energetic, these secondary electrons may excite or ionize other atoms near which they pass, giving rise to cascade of ionizing events. The damage to the cell is either direct or indirect ionization of the atoms which make up the DNA chain If there is molecular oxygen nearby the damage to the DNA-strand may be fixed. This phase is illustrated in Figure 11.

Figure 11: Direct and indirect action of charged particles. Image from [12]

Secondly there is a chemical phase. This is the period when these damage atoms and molecules interact with other cellular components in rapid chemical reactions. Ionization and excitation lead to the breakage of chemical bonds and formation of broken molecules known as free radicals. These are very reactive and start a series of reaction which damage the DNA.

Last there is a biological phase, which includes all subsequent processes. These begin with enzymatic reactions that acts on the residual chemical damage. The vast majority of lesions in for instance DNA are successfully repaired. Non-cancerous cells are able to reproduce even with slightly damaged DNA. The excited state of reproduction that cancerous cells are in means that a small amount of DNA-damage renders them incapable of reproducing. Some rare lesions fail to repair and these eventually lead to cell death (or mutation), a phenomenon used in radiotherapy. A cell takes time to die. After being irradiated they may undergo a number of mitotic diversions before dying. The radiation may also give rise to cell proliferation (increase in number of cells) and a very late effect is the appearance of a second tumour (called radiation carconogenesis). In the brain the injury in the central nervous system develop a few months to several years after therapy[12].

(18)

Figure 12: The area between the two curves are called the therapeutic window.

Another factor to take into considerations is the volume effect of the organ, or the seriality i.e. how the functional subunits is arranged. This effect is already mentioned in section 2.2.1. In the context of seriality, an organ is usually described as a series and parallel circuits. For instance the optic nerve is an organ with a high seriality, which would fit into the description of a series circuit. If tissue is destroyed so that the nerve is unconnected, the patient is likely to lose its sight. The lung is an organ with low seriality. If a part of it is destroyed there is still a probability that the function is maintained.

Figure 13: Figure illustrating the seriality. If a link is broken the first chain is no longer connected, but the second one is. Image from [23]

2.3.2 LQ-model

(19)

Figure 14: Explanatory LQ-plots. Image from [12]

According to the LQ-model, the surviving fraction of target cells after a dose per fraction d, (SFd), is

given as.

SFd= exp(−αd − βd2) (8)

where α and β is parameters fitted into the graph. An example of one such graph is in the left picture in Figure ??. One may compare the α_β-ratio between tissues. A high α_β-ratio implies that the tissue is early responding to radiation, while a low α_β-ratio implies that the radiation effects comes after some time. In the right picture in Figure ?? a comparison between the LQ-curves of a structure with early responding tissue and late responding tissue is made. One limitation of the LQ-model is that it works well for doses employed by conventional fractionated radiotherapy. For higher doses per fraction, as in stereotactic radiosurgery, the LQ-model is less accurate.

2.3.3 TCP and NTCP

Two common radiobiological models are TCP (Tumor Control Probability ) and NTCP (Normal-Tissue Complication Probability). TCP and NTCP describe probability that the cells of the tumour or normal tissue dies under the radiation treatment dose. Brahme [2] came up with the Poisson-based TCP-function.

T CP ( ¯d) = mT Y i=1 Pi where Pi= exp − Noexp(−(αdi+ βd2i) (9) Here N0 is the number of chlonogenic cells (a simplification is made assuming that the number of

chlonogens in each voxel is the same), α and β is derived from the LQ-model in Section 2.3.2, mT is the

number of voxels in the tumour, di is the dose in the i:th voxel and ¯d is the vector of all the di’s. The

TCP function is strictly concave and thus convenient to use in optimization [5]. The NTCP function is in K¨allman [18] described by N T CP ( ¯d) =1 − mR Y i=1 (1 − P_is)1/mR 1/s where Pi= exp − Noexp(−(αdi+ βd2i) (10) where di is the dose in voxel i, s ∈ (0, 1] is the relative seriality of the tissue , mR is the number of

voxels in the organ at risk and N0, α and β have the same function as in the TCP case. One issue with

(20)

2.3 Radiobiology 2 BACKGROUND 2.3.4 Equivalent Uniform Dose

Niemerko [21] made an assumption that two dose distributions are equivalent if they cause the same probability for tumor control, as demonstrated in Figure 15. Using this assumption and the TCP model as a basis he developed the Equivalent Uniform Dose-model, abbreviated as the EUD and defined in Equation 11. It reduces the complex three-dimensional dose distribution into a singe number. The EUD is essentially a norm [5] which takes the radiationsensitivity of the tissue into account.

EU D( ¯d) = −1 aln   1 |V | |V | X i=1 e−adi   (11)

Here, V is the set of all voxels in the structure of interest, a is a parameter specified for type of structure depending on the radiation response of the tissue (a ≤ 0 for risk structures and a < 0 for target structures) and di is the dose in the i:th voxel. Originally the EUD dealt with the effect of radiation in

tumour structure only. This was later simplified and extended to deal with organs at risk as well, by developing the generalized Equivalent Uniform Dose, the gEUD [22].

gEU D( ¯d) =   1 |V | |V | X i=1 dai   1 a (12)

The gEU D is a convex function for a ≥ 1 (concave for a ≤ 1) [5] and is thus more efficient to work with in optimization problems. While TCP and NTCP have many uncertain parameters, the EUD only has one parameter (although still uncertain). The EUD is based on both physical (dose) and biological criterias.

Figure 15: The EUD simplifies the dose distribution into one value. All distributions giving the same EUD are equivalent. Image from [23]

2.3.5 Hypoxia

(21)

Figure 16: Acute hypoxic, chronic hypoxic and oxic regions. Image from [12]

Oxygen makes the tumours more sensitive to radiation, i.e. it’s a strong radiosensitizer. It increases the effectiveness of a given dose of radiation by forming DNA-damaging free radicals. Tumor cells in a hypoxic environment may be as much as 2 to 3 times more resistant to radiation damage than those in a normal oxygen environment [12]. In Figure 17 one can see the increases resistance of the cells in a Hypoxic environment using the LQ-model.

Figure 17: Survival curve for cultured mammalian cells exposed to x-rays under oxic and hypoxic con-ditions. Image from [12]

To incorporate the resistance in the hypoxic cells one can modify the α and β parameter in the LQ-model and take into account the distance to molecular oxygen and the tissue sensitivity [27]. By calculating the OER as in Equation 13 one may assess the effect molecular oxygen has on the irradiated tissue, i.e. the enhancement of radiation damage.

Oxygen Enhancement Ratio = Radiation dose in hypoxia

Radiation dose in air (13) Based on the OER, the dose modification factor may be calculated.

f (r) = OERmax(k + pO2(r)) k + OERmax× pO2(r)

(14) Here OERmax is the maximum effect achieved in the absence of oxygen, k is the reaction constant and

pO2is the local oxygen tension, which depends on the distance r to oxygen. Oxygen effect only occurs if

(22)

2.4 Comparison between DVH, EUD, TCP and NTCP 2 BACKGROUND

α(r) = αoxic

f (r), β(r) = βoxic

(f (r))2 (15)

The modified TCP then becomes as in Equation 16.

T CP ( ¯d) = mT Y i=1 Pi where Pi= exp − N0exp(−(α(r)di+ β(r)d2i) (16) The only difference to the TCP stated in Equation 9 is the dependence of r in the α and β parameters.

2.4 Comparison between DVH, EUD, TCP and NTCP

(23)

3 OPTIMIZATION OF THE BEAM-ON TIME

3 Optimization of the Beam-on time

3.1 Problem set-up

The problem of finding the optimal treatment plan was divided into three parts. 1. Find the isocenter positions in the tumour.

2. Perform a dose calculation in order to find the dose rates in each voxel for a specific sector in a specific state.

3. Optimize the Beam-on time.

To find the isocenter position the grassfire-algorithm is applied and, when the isocenters are found, a simplified version of the TMR10 algorithm is used for the dose calculation. The grass-fire algorithm starts at the surface of the volume and goes inward to some sort of ”center”, defined by being most far away from the surface. When this center is found, the coordinates becomes the first isocenter. The volume covered by the isocenter is subtracted and the procedure is iterated. This is done until almost the entire volume is covered. Since it may be described as putting fire to the verge of a lawn, it’s named the grassfire-algoritm. Some other simplifications were that the number of isocenters and the isocenter positions were fixed throughout the optimization process. Also only the three of the different shapes of the isocenters were used in the procedure. All of them were spheres with a radius of either 4 mm, 8 mm or 16 mm. Fixating the positions in the isocenters renders a convex problem assuming the cost function is convex. All this was done to make the problem convex and decrease the number of degrees of freedom. The organs at risk and tumours were also divided into voxels to speed up the computations.

To calculate the dose in voxel i, di the sum below is used. This is a discretized version of the Fredholm

integral [6] and is similar the one used in Section 2.2.

di(¯ω) = J X j=1 K X k=1 L X l=1 ωjkldi,jkl (17)

Here ¯ω is the vector of all the beam-on times, j is an index over the isocenterpositions, k is a index over the sectors in the Gamma Knife (in total 8 sectors) and l is a index over the three collimator sizes. di,jkl

is the doserate in voxel i from sector k using collimatorsize l when the isocenter is in position j and ωjkl

is the beam-on time of sector k using collimatorsize l when the isocenter is in position j. Using this estimation of the dose the treatment plan will be found by optimizing over ωjkl. The problem will have

J KL degrees of freedom where J = 3, K = 8 and L is the total number of isocenters. Throughout this thesis, the sum in Equation 17 wil be written as di(¯ω) =Pµωµdi,µ, where µ = 1, 2, 3, . . . J KL.

3.2 The Radiotherapeutical model-based on DVH

A DVH-based objective function is hard to solve due to the non-convexity of the DVH. Instead a sim-plification of the description of the DVH is made by using a nonlinear penalty function penalizing when the dose does not fulfill the constraints. Three such constraints are implemented.

3.2.1 Maximum and minimum dose-constraints

To penalize when the dose exceeds the maximum dose, dmax, the function P m

i=1(di(¯ω) − dmax)2+ is

implemented in the objective, where (x)+ = max(x, 0), m is the number of voxels in the structure

and di is calculated as in Equation 17. The objective is to push the dose down to dmaxusing a squared

penalization factor over all the voxels in the structure. The minimum dose constraint is treated similarly.

m

X

i=1

(dmin− di(¯ω))2+

(24)

3.3 The Radiobiological model 3 OPTIMIZATION OF THE BEAM-ON TIME 3.2.2 Dose-Volume constraints

The Dose-Volume constraints, or DV-constraints, are specified as V ({v : dv≥ dr}) < Vr, or the volume

receiving the dose drshould be less than Vr. To implement this one needs to find dpsuch that V (dp) = Vr.

Hence we have first the required dose dr and the present dose dp at Vr. This is illustrated in Figure 18.

By penalizing all doses between dr and dp as in Equation 18 one may pull the curve closer to dr. This

method was first suggested by Bortfeld [1].

H(di(¯ω) − dr)H(dp− di(¯ω))(di(¯ω) − dr)2 (18)

where H(·) is the Heaviside function. This is done similarly for each DV-constraint. This type of constraints are at the moment the most clinically used one.

Figure 18: Illustration of how dp is found.

3.2.3 The resulting DVH model

The resulting and somewhat complicated model is then:

(N LP )                  min P t∈T λt mt Pmt i=1 (dmin− di(¯ω))2++ P p∈DV CH(di(¯ω) − dp)H(dr− di(¯ω))(di(¯ω) − dp)2+ (di(¯ω) − dmax)2+ +P r∈R λr mr Pmr i=1 (di(¯ω) − dmax)2++ P p∈DV CH(di(¯ω) − dp)H(dr− di(¯ω))(di(¯ω) − dp)2 + λP µωµ s.t. di(¯ω) =Pµωµdi,µ, i = 1, 2, 3, . . . ωµ≥ 0 (19) where the λ’s is the weight on the treatment time, organ or tumour, H(·) is the Heaviside function T are the tumours, R are the riskorgans and (x)+= max(x, 0). In short, a minimum dose and a DV-constraint

is implemented for the tumours and a maximum dose and DV-constraint is implemented for each organ at risk. Also, the total treatment time is penalized by the last term in the objective function. There is also a possibility to specify the importance of each organ by using the λ’s to weight the organs between one another. If a part of the tumour is hypoxic, this are is treated as a separate tumour with other constraints. The square term is chosen mostly because the full objective function remains convex[28]. One may also penalize the constraints harder by changing the quadratic penalization function into an exponential function, (·)4or another convex function.

3.3 The Radiobiological model

3.3.1 gEUD-based objective

(25)

3.4 GUI 3 OPTIMIZATION OF THE BEAM-ON TIME (N LP )          min P t∈T λt(gEU D0− gEU D(¯ω))2++ P r∈Rλr(gEU D(¯ω) − gEU D0)2++ λ P µωµ s.t. gEU D(¯ω) = (P|V | i=1di(¯ω)a) 1 a di=P_µωµdi,µ, i = 1, 2, 3, . . . ωµ≥ 0 (20)

Since the gEU D isn’t very intuitive, gEU D0 is derived from the desired DVH-curve. The value of the

parameter a is uncertain, hence the model is solved and tested for the a parameter in an interval of 20% around the mean to study the robustness of the optimized plan.

3.3.2 TCP/NTCP-based objective

The TCP has a great advantage of being concave, which is not the case for the NTCP function. Although by applying the transformation h(z) = −ln(1 − z) yields a convex NTCP [15], [24].

(N LP )          min P t∈T λt T CP0(T CP0− T CP (¯ω)) 2 ++ P r∈R λr ln(1−N T CPo)(−ln(1 − N T CP (¯ω)) + ln(1 − N T CP0)) 2 + +_TλP µωµ s.t. di(¯ω) =Pµωµdi,µ, i = 1, 2, 3, . . . ωµ≥ 0 (21) Here T is a normalizing constant, which is the time it takes to irradiate the tumour with just the 8mm collimator, serving the purpose of making the cost function dimensionless. Also here, the parameters are uncertain so a sensitivity analysis is performed for a 20% interval around the mean. The model is very sensitive to the initial solution. In order to give an initial estimate, an absolute value of a plot is chosen using a simple solution to the DVH-model as an initial estimate.

3.4 GUI

(26)

3.4 GUI 3 OPTIMIZATION OF THE BEAM-ON TIME To be able to optimize in the gui some initial steps have to be taken. The first one is made in the ”Con-sidered organs” panel. By clicking the button ”New Patient” a window pops up containing the names of all the outlined structures in the patient and an individual number to each structure. By writing a vector in the editbox containing the numbers to the target structures and clicking ”Ready”, the program starts to take the info from the DICOM-structures and create all the info needed to continue with the initial steps and the optimization and saving the info as s.mat in the patient folder. Hence, the program can only upload the file s.mat the next time one wishes to investigate the patient, to gain speed. When this is done, one needs to define organs at risk and which target to consider. This is done by clicking the ”New System” button. A plot is displayed with first all the targets (to chose targets), then all the outlined organs (to chose organs at risk). Then if there is a large problem with many tumours one can split them into smaller subsytems by clicking the target and organs at risk and pressing ”Next Sytem”. This division will later affect the dose-matrix by making it sparse to speed up calculations during the optimization. If one wishes to add a hypoxic area in the tumour, one can outline it by pressing ”Add hypoxic area”.

In the isocenter-panel the isocenters are defined. First the Cut-off distance is specified. The Cut-off distance is the smallest distance to another isocenter and in this case a sort of tolerance. When pressing ”Pack spheres” the grassfire-algorithm starts to calculate the isocenters. To gain degrees of freedom only the 4mm spheres are used. The number of isocenters will be displayed in the black box and they will be plotted in the tumour in the graph.

Options for optimization, such as maximum number of iterations and function evaluations along with the tolerance, can then be specified. To speed up the calculations decreasing tolerances may be specified in a vector. MATLAB then solves the optimization for each tolerance using the results from the previous as an inittial solution to the next. This is in many cases preferred to make computations faster. By pressing the ”Set parameters” button a window will pop-up with a button for each target and organ at risk in the problem. There is also a box for specifying the weight on the time and a button for showing some parameters for different organs in the brain. By pressing the button for each organ parameters will show up. One may also specify model and algorithm in the pop-up boxes. The model ”TCP hypoxia” is a bit different from the others. Here one can instead outline blood vessels and the parameters α and β is then perturbed according to Equation 15. The p0(r) is, for simplicity assumed to be a linear function, where p0(0) is the oxygen tension in blood (p0=100mmHg), and r such that p0(r) = 0 is the diffusion distance

of blood in a muscle (rd = 39µm). The equation approximating the oxygen tension is displayed below.

Other parameters used when perturbing α and β is OERmax and k, which is set to 3 and 2.5mmHg

respectively. p0(r) = ( p0(1 −|r|_r_d) if |r| ≤ rd 0 if |r| > rd (22) Then every initial step is made and one can proceed to optimization by clicking ”Start Optimization”. Then dose-calculation is initiated. When this is done a simplified system is constructed containing only targets and organs at risk and the simplified dose-matrix. The assumption is made that a irradiat-ing a target in a subsystem doesn’t affect the dose in other subsystems considerably. The dose-matrix (diµ)i=1,..,m,µ=1,...,J KL used when calculating the dose in the voxels (in Equation 17) is made sparse

by taking the mean over a specified amount of rows in the columns describing the doserate from the isocenters in the other systems. This makes calculations faster during the optimization. Every step in the optimization is displayed in the MATLAB command window.

(27)

3.4 GUI 3 OPTIMIZATION OF THE BEAM-ON TIME

Figure 20: The interface of the GUI for isodoseareas.

The N0-parameter may be estimated by first enter the values for α and β in the target and then pressing

”DVH using LGK”. N0 is then calculated assuming the Gamma Knife treatment killed 90% of the

tumour. Also the EU D0may be estimated two approaches. First options is to set EU D0 to zero before

(28)

4 RESULTS AND DISCUSSION

4 Results and Discussion

Some test-cases were provided and two of them are presented here. These cases are patients already treated in the Gamma Knife. The plans were derived using the interior-point algorithm and all three models. The DVH and EUD-based models were also solved in the case of adding a hypoxic region in the tumour, and the TCP-model, which takes hypoxia into account by perturbing the α and β-parameters where solved.

4.1 Case 1: One tumour

The first test-case is a tumour located close to the brainstem and the cochlea. The target and the organs at risk and the resulting DVH from the Gamma Knife is below. EU D0 is estimated using the results

from solving the DVH-based model.

Figure 21: Mesh of the head, where tumour and organs at risk are outlined and the DVH from the resulting Gamma Knife dose distribution.

(29)

4.1 Case 1: One tumour 4 RESULTS AND DISCUSSION Parameter Tumour Brainstem Cochlea

λ 1 0.5 0.1 dmin 14 - -dmax 28 15 10 DVC (18, 0.7), (23, 0.3) (8, 0.2) -a -8 7 5 EU D0 18.01 8.772 8.351 α 0.22 0.0491 0.0878 β 0.0272 0.0234 0.0293 s - 1 0.64 N0 10 5 5 T CP0 100% - -N T CP0 - 1% 1%

The Cut-off distance 1.5 mm generated 14 isocenters. Below is a plot of the chosen hypoxic area. This hypoxic area has the same constraints as the tumour only the value of dmin is set to 25Gy. The tumour

size is 4.9cm3 _{and the hypoxic region is placed at the top of the tumour. The hypoxic region is usually}

in the center of the target, but is placed in the top here to display more clearly if there will be any difference in the dose in the lower and upper part of the tumour.

Figure 22: The hypoxic region of the tumour

4.1.1 Results DVH without hypoxic region

In this case λtimeis set to 10−4. Solving the optimization problem yields a value of first order optimality

(30)

4.1 Case 1: One tumour 4 RESULTS AND DISCUSSION

Figure 23: DVH of the target and organs at risk where each constraints is plotted.

Above is a plot of the DVH. The thicker lines are the results from the DVH model and the thinner lines is the results from Leksell Gamma Plan R_{. The horizontal lines and dots are the constraints set. The} result follow the defined soft constraints as one can see in the DVH above and there is a coverage of 100%. Also the gradient of the dose in the normal tissue is steep and organs at risk are steep. The dose in the normal tissue is much higher than en the Leksell Gamma Plan R _{plan due the the fact that it is} not penalized.

Figure 24: Histogram of the Beam-on times.

(31)

Figure 25: Isodose of the resulting values.

In the isodoseplot the contours of the dose above 14 and 20 Gy is plotted along with contours of the organs and for the target. One can clearly see that most of the dose is distributed so that the brainstem and the cochlea are minimally affected. One can also see this behavior in the isosurfaceplot in Figure 26.

Figure 26: Isodosesurface of the resulting weights.

Testing robustness in λ

(32)

Figure 27: The DVH for the different values of the parameter λT arget. Here the DVH from the previously

solved plan is the thin line.

One can see that the model is fairly robust when varying λT arget.

Testing the algorithms

To test the other algorithms the DVH-problem was solved for the three other algorithms as well.

(33)

Figure 29: Results using the SQP algorithm. The resulting function value was 1.08 × 10−2

Figure 30: Results using the trust-region reflective algorithm. The resulting function value was 6.12×10−3 All algorithms give similar results. By comparing the function value (SQP gave a function value of 1.08 × 10−2) the trust-region reflective provided the best result, although this algorithm is very slow (total iteration-time was 3366 seconds compared to the other algorithms which solved the problem in around a couple of minutes). The SQP and the interior-point algorithm also have the advantage that it can recover from infinite and not defined values of the objective function, which is needed in the TCP-model.

(34)

4.1 Case 1: One tumour 4 RESULTS AND DISCUSSION calculation took 287 seconds and the beam-on time was 114 minutes.

Figure 31: DVH-curve of the target and organs at risk where each constraints is plotted. The coverage was 100% of the target and the fall-off dose of the gradients are steep. One can also clearly see the small bump in the DVH for the tumour due to the influence from the hypoxic region. The coverage of the hypoxic region is 87%. By penalizing this higher one may achieve a better coverage, but there will be a trade-off between the minimum dose in the hypoxic region and of the maximum dose in the tumour and organs at risk as well as the normal tissue. The integral dose in the target is 205mJ, compared to the non-hypoxic case when it was 186mJ.

Figure 32: Isodose of the result. Comparing the the non-hypoxic case the dose is higher in the upper part of the tumour.

(35)

Figure 33: Isodosesurfaces for the dose above 14 and 10 Gy.

One remarkable observation here is the increased beam-on time. Setting the weights for all the organs at risk to 0 and redo the calculations yields a beam-on time of 66 minutes and the DVH below.

Figure 34: DVH of the resulting Beam-on times in the case when no weight is put on the organs at risk. Due to the constraints on the organs at risk, less contribution will be given from the sectors, whose beams pass through the organs at risk. Due to the increased beam-on time when giving higher dose to a part of the tumour which is just in the vicinity of these organs at risk a lot less contribution will be given, hence the beam-on time increases drastically.

4.1.3 Results EUD without hypoxic region

(36)

4.1 Case 1: One tumour 4 RESULTS AND DISCUSSION Organ EUD

Tumour 18.01 Cochlea 8.76 Brainstem 8.35

Figure 35: DVH of the resulting weights.

Here the resulting dose is somewhat lower then in the previous case due to he fact that either the maximum nor minimum dose is controlled. Even though minimum and maximum doses aren’t taken into consideration in this case, the plan fulfills the DVH-constraints reasonably well.

Figure 36: Isodoses of the target and organs at risk

(37)

Figure 37: Isodosesurfaces of the target and organs at risk.

Testing robustness in a

The parameter a was varied in an interval of 20% around the nominal value and solved in the optimization for the target and brainstem to test the robustness. The basic case was when all parameters were set to the nominal value of the parameter, i.e. the case previously solved.

(38)

4.1 Case 1: One tumour 4 RESULTS AND DISCUSSION 4.1.4 Results EUD with hypoxic region

Due to the hypoxic region being a part of the tumour and they will be regarded as two seperate organs the EU D0 will be modified. The EU D0 for the tumour and hypoxic area was set to 19.19 and 26.03

respectively. The number of iterations was 1145 and the first order optimality of 1.47 × 10−2. The calculation time was 142 seconds and the beam-on time was 50 minutes. The resulting EUD is displayed in the table below.

Organ EUD Tumour 19.18 Hypoxic area 26.02 Cochlea 8.77 Brainstem 8.35

One can see here that no resulting EUD deviates much from the target EUD.

(39)

Figure 40: Isodoses of the target and organs at risk.

Here one may notice a higher dose in the upper part of the tumour when comparing subplot one and two with the non-hypoxic case, which is the hypoxic region.

Figure 41: Isodosesurfaces of the target and organs at risk.

4.1.5 Results TCP without hypoxic region

Here λT umourwas set to 100. The number of iteration needed to solve the problem was 2631 and it took

317 seconds. The value of first order optimality was 1.786 × 10−1 and the total beam-on time was 492 minutes. The resulting value of TCP and NTCP is shown in the table below.

Organ TCP Tumour 99%

(40)

The coverage was at 95%. Even thought the dose in the organs at risk is less then in the previous cases the dose in the tumour is very high (even too high). The same phenomenon is mentioned in [16]. Since the TCP don’t have any constraints on the maximum dose at all, the dose can grow very high. It will change if the healthy tissue is taken into account in the optimization.

Figure 43: Isodosearea of the result were doses above 14 and 10 Gy are marked.

(41)

Figure 44: Isodosesurfaces of the result were doses above 14 and 10 Gy are marked.

Including normal tissue in optimization

The theory about including the healthy tissue to avoid a too high dose in the target is tested. The parameters for the healthy tissue is set to α = 0.0499, β = 0.0238, s = 0.64 and λ = 0.1. An increase was made in λtimeto 0.05. Below is the DVH of the resulting plan.

Figure 45: DVH when including the normal tissue.

(42)

Figure 46: Isodoses of the results.

The table below displays the resulting TCP and NTCP. The NTCP in the healthy tissue and in the organs at risk is very high due to some parameters having erroneous values.

Organ TCP Tumour 76% Organ NTCP Cochlea 32% Brainstem 28% Normal tissue 80%

Testing robustness in α, β and N0

Also here a robustness check were performed giving the following results. α and β were varied using three values, one for 0.8 and one for 1.2 of the value and one using the expected value. For N0, N0,T umour

varied between 100 and 1000 number of chlonogenes keeping N0 for organs at risk at 5, and N0,Cochlea

and N0,Brainstem varied between 10 and 100 number of chlonogens keeping N0,T umour at 1000 number

(43)

Figure 47: The DVH for the different values of the parameter α. The thin line is the original result.

(44)

Figure 49: The DVH for the different values of the parameter N0. The thin line is the original result.

4.1.6 Results TCP with hypoxic region

(45)

To solve the problem 1815 iteration were performed and took 498 seconds. The first order optimality was 0.52 × 10−6 and the total beam-on time was 259 minutes. The λ’s was same as in the TCP-problem without hypoxic region.

Organ TCP Tumour 99%

Organ NTCP Cochlea 2.8% Brainstem 5.3%

(46)

Figure 52: Isodose of the result were doses above 14 and 20 Gy are marked.

The isodose plots also show no bigger difference between this case an the non-hypoxic case. Also here they may be deceptive.

Figure 53: Isodosesurfaces of the result were doses above 14 and 10 Gy are marked.

4.1.7 Comments

(47)

4.2 Case 2: Two tumours 4 RESULTS AND DISCUSSION

Figure 54: A TCP and NTCP graph.

However, when an initial estimate can return a value of NTCP and TCP both much larger than zero (as in the figure below), the optimization will work as it should. This is also just a numerical problem due to the theoretical fact that these curves are not likely to reflect the true clinical situation. The tumour is in general more sensitive to irradiation than the healthy tissue and it’s always the tumour which is irradiated with the highest dose.

4.2 Case 2: Two tumours

(48)

Figure 55: Mesh of the head, where tumour and organs at risk are outlined and the DVH from the resulting Gamma Knife dose distribution.

Parameter Right Frontal Me Left Parietal Po Right Thamalic F Foot

λ 1 1 0.5 1

dmin 15 15 -

-dmax 30 30 15 15

DV C (20 0.3),(17 0.7) (20 0.3),(17 0.7) (10, 0.5) (10, 0.5)

The estimate if N0yields a value below 1. Instead in these calculations N0 is set to 10 for the tumours

and 5 for the organs at risk. Here EU D0is calculated as in Case 1, and λtimeis set to 10−4. 33 isocenter

is distributed using the grassfire-algorithm with cut-off distance of 1.7mm.

4.2.1 Results DVH without hypoxic area

λT ime= 5 × 10−5. The optimization took 1051 seconds and 1709 iterations. The value of the first order

(49)

Figure 56: DVH of the resulting values. The thin line is the result from the clinical plan. The DVH show a steep fall-off dose and a coverage of 100%. The brainstem get a slightly higher dose than the maximum constraint, although by penalizing it more will yield a better fit.

Figure 57: Isodoses of the dose.

(50)

Figure 58: Isodosesurface of the resulting dose.

4.2.2 Comments

(51)

5 CONCLUSION

5 Conclusion

The DVH-based model worked well according to plan and the solution is good when comparing to the clinical plan and regarding the defined soft constraints. The EUD-based model did also work well as well and followed the soft constraints reasonably. Both models were easy to work with by giving harder constraints or weighting the organs. There were no problems in implementing these models and they all performed well using any MATLAB-solver algorithm. The main concern was instead when the beam-on time was minimized. This was implemented as a L1-problem, so that it could keep the beam-on times low and reduce the incidence of two or more collimators irradiating from one sector. However, by penal-izing the beam-on time too much may force many sectors to be turned off, which give plans with poor selectivity. Hence one must find a balance in the time-penalty.

The TCP/NTCP-based model however was a bit more troublesome implementing. Since the gradient of TCP and NTCP are larger than zero only in a small region, choosing a initial point in this region is vital for the numeric optimization. If the initial point is too low, the TCP will have a low gradient and will be iterated down to zero. If on the opposite the initial point is too high the initial estimate will yield an infinite objective. Also one must weight the tumours much higher than the organs at risk and the time due to the fact that the TCP belongs to a compact interval (TCP∈ [0, 1]). Even though the relative difference to the reference value is taken, the convex transform of the NTCP is by a logarithmic function, causing small differences in ω to become large when calculating −ln(1 − N T CP ). Hence the problem may once again iterate down to ¯ω = ¯0 if the weight on the TCP is too low.

The method of guiding a higher dose to hypoxic regions was proven to be successful in the case of DVH and EUD. Also the method in the TCP-case were proven to be successful.

5.1 Future work

(52)

A OPTIMIZATION ALGORITHMS

A

Optimization algorithms

Matlabs fmincon-function was used to solve the models. The algorithms used in fmincon can either be the Active-set algorithm, Sequential Quadratic Programming (SQP), Interior-point method or the Trust-region relective method. A short description of these algorithms will be given in this section and the derivations are acquired from [13]. Suppose the problem to be solved is:

(N LP )      min f (x) s.t gi(x) ≤ 0, ∀i ∈ I x ∈ Rn (23)

The Lagrangian function is defines as L(x, λ) = f (x) − λTg(x). The first order optimality condition is ∇L(x, λ) = 0 and the second order optimality condition is that ∇2xxL(x) is positive definite. For a

solution to be optimal, both have to be fulfilled.

A.1 Active-set algorithm

The active-set-algorithmis is a fast algorithm since it can take large steps. The algorithm is effective on some problems with nonsmooth constraints. It is not a large-scale algorithm and it generates feasible points in each step.

For x to fulfill the first order optimality condition (∇L(x, λ) = 0), Equation 24 have to be satisfied. ∇xL ∇λL =∇f (x) − A(x) T_λ g(x) =0 0 (24) This is also called the Karusch-Kuhn-Tucker conditions. Here A(x) = (∇g1(x), ∇g2(x), . . . ) and it is

assumed thar A(x) has full row-rank. This is solved using the Newton iteration, setting x+ = x + p, λ+= λ + ν in Equation 24 and make a second order Taylor expansion. The resulting equations is.

∇2 xxL(x, λ) −A(x)T A(x) 0 p ν =−∇f (x) + A(x) T_λ −g(x) (25) This can be simplified into Equation 26.

∇2 xxL(x, λ) A(x)T A(x) 0 p −λ+ =−∇f (x) −g(x) (26) Compare this to the first-order optimality condition for a Quadratic Problem solved by a Newton step in Equation 27. (QP ) ( min 1₂xT_{Hx + c}T_x s.t Ax = b ⇒ H AT A 0 p∗ −λ∗ = −Hx + c Ax − b (27)

By using these similarities one obtain a new optimization problem for the step p below

(QP )      min 1₂pT_∇2 xxL(x, λ)p + ∇f (x)Tp s.t A(x)p = −g(x) p ∈ Rn (28)

Now, the optimal step need to be found using the Newton step setting p = ¯p + p0. Assume that the feasible point ¯p is known. Guess that the constraints active at ¯p are also active at the optimal p∗. Let A = {I : aT

I(x)¯p = −gI(x)}, i.e. the constraints active at ¯p. Let W ⊆ A be such that AW(x) have full

row rank. Keep temporarily the constraints in W active, i.e. solve

(EQPW)      min 1₂(¯p − p0)T∇2 xxL(x, λ)(¯p − p0) + ∇f (x)T(¯p − p0) s.t A(x)p0 = 0 p ∈ Rn (29)

(53)

A.2 Sequential Quadratic Programming A OPTIMIZATION ALGORITHMS ∇2 xxL(x, λ) AW(x)T AW(x) 0 p∗ −λ∗ = −∇ 2 xxL(x, λ)¯p + ∇f (x) 0 (30) Since the (QP)-problem solved is equality constrained, no further description of how to include in-equalities is made. To summarize the Active-set algorithm one iteration is described. Given a fea-sible ¯x to Equation 23, a feasible ¯p to Equation 29 and W such that AW(¯x) has full row rank and

AW(¯x)¯p = −gW(¯x). 1. Solve∇ 2 xxL(x, λ) AW(x)T AW(x) 0 p∗ −λ∗ = −∇ 2 xxL(x, λ)¯p + ∇f (x) 0

2. If λW ≥ 0 then ¯p is optimal and the solution for the (QP)-problem is found. Let p = ¯p + p∗.

Otherwise new (QP)-iteration. 3. Let x = ¯x + p and solve Equation 26.

4. If λ+≥ 0, x is optimal. Otherwise, make a new iteration.

A.2 Sequential Quadratic Programming

SQP satisfies bounds at all iterations and can recover from NaN or Inf results, although it is not a large-scale algorithm. Assume the second order optimality condition is fulfilled. The proceedings are the same as in the Active-set method until Equation 29. Then instead of following the active constraints Assume A(x) has full row rank. One iteration with the SQP-solver for a nonlinear problem is then:

1. Complute optimal solution p and multiplier vector λ+ to

(QP )      min 1₂pT_∇2 xxL(x, λ)p + ∇f (x)Tp s.t A(x)p = −g(x) p ∈ Rn 2. Change x to x + p, and λ to λ+

The difference between SQP and active-set method is that in the subproblem for the active-set method you only consider the working sets, i.e. the set of constraints which is currently active, and then solve for the other constraints. While in the SQP subproblem you solve for all constraints.

A.3 Interior-point method

The Interior method has the advantage that it can be used on very large problems. However, the solutions can be slightly less accurate than those from other algorithms because of the barrier function, which keeps the solution away from constraint borders. Although for most practical purposes, this inaccuracy is usually quite small. The interior-point method satisfies bounds at all iterations, and can recover from NaN or Inf results. Also it can use special techniques for large-scale problems.

First define the Barrier function as in Equation 31.

Bµ(x) = f (x) − µ m

X

i=1

ln(gi(x)) where µ > 0 (31)

Bµ(x) keep the solution away from the constraint border. A necessary conditions for a minimizer of

Bµ(x) is fulfilling the first-order optimality condition, i.e. ∇Bµ(x) = 0 where

∇Bµ(x) = ∇f (x) − µ m X i=1 1 gi(x) ∇gi(x) = ∇f (x) − µA(x)TG(x)−1e (32)

(54)

A.4 Trust-region-reflective A OPTIMIZATION ALGORITHMS ∇f (x(µ)) − A(x(µ))T_{λ(µ) = 0} λ(µ) − _g(x(µ))µ = 0, i = 1, . . . , m⇐⇒ ∇f (x(µ)) − A(x(µ))T_{λ(µ) = 0} gi(x)λ(µ) − µ = 0, i = 1, . . . , m (33) Taking a Newton step and making a second order Taylor expansion yields Equation 34.

∇2 xxL(x, λ) A(x)T ΛA(x) −G(x) ∆x −∆λ =∇f (x) − A(x) T_λ G(x)λ − µe (34) where Ω equals to At last, given µ > 0, x such that g(x) > 0 and λ > 0, one iteration with the iteror-method is: 1. Compute ∆x, ∆λ from ∇2 xxL(x, λ) A(x)T ΛA(x) −G(x) ∆x −∆λ = −∇f (x) − A(x) T_λ G(x)λ − µe

2. Choose suitable steplength α such that g(x + α∆x) > 0 and λ + αλ > 0. 3. Change x to x + α∆x and λ to λ + α∆λ.

4. If (x, λ) sufficiently close to (x(µ), λ(µ)), reduce µ.

A.4 Trust-region-reflective

The main idea with the Trust-Region-reflective method is find a Trust-region at the current point and improve the function value [4]. If you have a point x in a space and a function f (x) you want to minimize, then in a small neighborhood N of this point x you can estimate f with a simpler function q(x) instead. This neighborhood is called a Trust-Region. By knowing the Trust-Region and the approximate function q(x) to f (x) one may instead find a trial step s by computing the optimal to the problem mins{q(s), s ∈ N }. The current point is then updated to x + s if f (x + s) < f (x), otherwise x is

unchanged, N is shrunk and the trial-step calculations is repeated. The question is now how to find q(x) and N . The Trust-Region subproblem is usually stated as (QP)-problem similar to Equation 29.

(QP ) (

min 1₂sTHs + sTg

s.t. kDsk ≤ ∆ (35)

where g is the gradient of f and H is the Hessian of f at the current point x, D is a diagonal scaling matrix, depending on the constraints, ∆ ∈ R+ _{and k · k is the 2-norm.}

An iteration with the Trust-Region-reflective algorithm is 1. State the Trust-Region subproblem

2. Solve (QP ) ( min 1 2s T_{Hs + s}T_g s.t. kDsk ≤ ∆ 3. If f (x + s) < f (x), then x = x + s. 4. Adjust ∆.

(55)

REFERENCES REFERENCES

References

[1] T. Bortfeld. Optimized planning using physical objectives and constraints. Seminars in Radiation Oncology, 9:20–34, 1999.

[2] A. Brahme. Optimized radiation therapy based on radiobiological objectives. Seminars in Radiation Oncology, 9:35–47, 1999.

[3] A. Brahme, J.E. Roos, and I. Lax. Solution of an integral equation encountered in radiation therapy. Physics in Medecine and Biology, 27:1221–1229, 1982.

[4] H. Byrd, R. Schnabel, and G. Shultz. A trust region reflective algorithm for nonlinearly constrained optimization. SIAM Journal of Numerical Analysis, 24:1152–1170, 1987.

[5] B. Choi and J.O. Deasy. The generalized equivalent uniform dose function as a basis for intensity-modulated treatment planning. Physics in Medecine and Biology, 47:3579–3589, 2002.

[6] A. Chvetsov, J. Dempsey, and J. Palta. Optimization of equivalent uniform dose usinf the l-curve criterion. Physics in Medecine and Biology, 52:5973–5984, 2007.

[7] J.O. Deasy. Multiple local minima in radiotherapy optimization problems with dose-volume con-straints. Medical Physics, 24:1157–1161, 1997.

[8] Internal document. The convolution algorithm in leksell gammaplan 10. WhitePaper Article no: 018881.02, September 2011.

[9] Internal document. Inverse planning in leksell gammaplan 10. White Paper Article no: 018880.02, September 2011.

[10] Internal document. A new tmr dose algorithm in leksell gammaplan. White Paper Article no: 1021357.00, July 2011.

[11] A. Donnet, D. Valade, and J R´egis. Gamma knife treatment for refractory cluster headache: prospec-tive open trial. Journal of Neurology, Neurosurgery and Psychiatry, 76:218–221, 2005.

[12] G. Gordon Steel et al. Basic Clinical Radiobiology. Arnold, 2002.

[13] I Griva, S. Nash, and A. Sofer. Linear and Nonlinear Optimization. SIAM, Philadelphia, 2009. [14] J.M. Herman, J.H. Petit, P. Amin, Y. Kwok, P.R. Dutta, and L.S. Chin. Repeat gamma knife

radiosurgery for refractory or recurrent trigeminal neuralgia: treatment outcomes and quality-of-life assessment. International Journal or Radiation Oncology Biology and Physics, 59:112–1226, 2004. [15] A. Hoffmann, D. den Hertog, A. Siem, J. Kaanders, and H. Huizenga. Convex reformulation of

biologically-based multi-criteria intensity-modulated radiation therapy optimization including frac-tionation effects. Physics in Medecine and Biology, 53:6345–6362, 2008.

[16] L. Jones and P. Hoban. A comparison of physically and radiobiologically based optimization for imrt. Medical Physics, 29:1447, 2002.

[17] Y. Kwon and C.J. Whang. Stereotactic gamma knife radiosurgery for the treatment of dystonia. Stereotactic and Functional Neurosurgery, 64:222–227, 1995.

[18] P. K¨allman, A. ˚Agren, and A. Brahme. Tumor and normal tissue response tofractionated non uniform dose delivery. International Journal of Radiation Biology, 62:249–262, 1992.

[19] L. Leksell. The stereotaxic method and radiosurgery of the brain. Acta Chirurgica Scandinavica, 102:316, 1951.

(56)

REFERENCES REFERENCES [22] A. Niemerko. A generalized concept of equivalent uniform dose. Medical Physics, 26:1100, 1999. [23] A. Niemierko. Image-Guided IMRT, chapter Biological optimization. Springer, Berlin Heidelberg,

2006.

[24] E. Romeijn, J. Dempsey, and J. Li. A unifying framework for multi-criteria fluence map optimization models. Physics in Medecine and Biology, 49:1991–2013, 2004.

[25] J. R´egis, F. Bartolomei, M. Hayashi, and P. Chauvel. What role for radiosurgery in mesial tomporal lobe epilepsy. Zentralblatt fuer neurochirurgie., 63:101–105, 2002.

[26] I Toma-Dasu, A. Dasu, and J. Fowler. Should single or distributed parameters be used to explain the steepness of tumour control probabolity? Physics in Medecine and Biology, 48:387–397, 2003. [27] I Toma-Dasu, L. Uhrdin, J.and Antonovic, A. Dasu, A. Nuyts, P. Dirix, K. Haustermans, and

A. Brahme. Dose prescription and treatment planning based on fmiso-pet hypoxia. Acta Oncologica, 51:222–230, 2012.

[28] S. Webb. The physical basis of imrt and inverse planning. Brittish journal of radiobiology, 76:678– 689, 2003.

(57)

(58)

TRITA-MAT-E 2014:02 ISRN-KTH/MAT/E—14/02-SE

Development of a test bench for Gamma Knife Optimization

Development of a test bench for

Gamma Knife Optimization

V I C T O R I A S V E D B E R G

Development of a test bench for

Gamma Knife Optimization

V I C T O R I A S V E D B E R G

Contents

1

Introduction

2

Background

2.1

Radiosurgery

2.2

Radiotherapy

2.3

Radiobiology

2.4

Comparison between DVH, EUD, TCP and NTCP

3

Optimization of the Beam-on time

3.1

Problem set-up

3.2

The Radiotherapeutical model-based on DVH

3.3

The Radiobiological model

3.4

GUI

4

Results and Discussion

4.1

Case 1: One tumour

4.2

Case 2: Two tumours

5

Conclusion

5.1

Future work

A

Optimization algorithms

A.1

Active-set algorithm

A.2

Sequential Quadratic Programming

A.3

Interior-point method

A.4

Trust-region-reflective

References