Institute of Physics and Engineering in Medicine Phys. Med. Biol. 60 (2015) 2179–2194

Physics in Medicine & Biology doi:10.1088/0031-9155/60/6/2179

A new optimization method using a compressed sensing inspired solver for real-time LDR-brachytherapy treatment planning C Guthier1, K P Aschenbrenner1, D Buergy2, M Ehmann2, F Wenz2 and J W Hesser3 1

  Department of Experimental Radiation Oncology, Medical Faculty of Mannheim, Heidelberg University, 68167 Mannheim, Germany 2   Department of Radiation Oncology, Medical Faculty Mannheim, Heidelberg University, University Medical Centre Mannheim, 68167 Mannheim, Germany 3   Department of Experimental Radiation Oncology, Medical Faculty Mannheim, Heidelberg University, 68167 Mannheim, Germany and IWR, Heidelberg University, 69126 Heidelberg, Germany E-mail: [email protected], Katharina.Aschenbrenner@ medma.uni-heidelberg.de, [email protected], [email protected], [email protected] and [email protected] Received 11 June 2014, revised 29 December 2014 Accepted for publication 6 January 2015 Published 16 February 2015 Abstract

This work discusses a novel strategy for inverse planning in low dose rate brachytherapy. It applies the idea of compressed sensing to the problem of inverse treatment planning and a new solver for this formulation is developed. An inverse planning algorithm was developed incorporating brachytherapy dose calculation methods as recommended by AAPM TG-43. For optimization of the functional a new variant of a matching pursuit type solver is presented. The results are compared with current state-of-the-art inverse treatment planning algorithms by means of real prostate cancer patient data. The novel strategy outperforms the best state-of-the-art methods in speed, while achieving comparable quality. It is able to find solutions with comparable values for the objective function and it achieves these results within a few microseconds, being up to 542 times faster than competing state-of-the-art strategies, allowing real-time treatment planning. The sparse solution of inverse brachytherapy planning achieved with methods from compressed sensing is a new paradigm for optimization in medical physics. Through the sparsity of required needles and seeds identified by this method, the cost of intervention may be reduced. 0031-9155/15/062179+16$33.00  © 2015 Institute of Physics and Engineering in Medicine  Printed in the UK

2179

C Guthier et al

Phys. Med. Biol. 60 (2015) 2179

Keywords: brachytherapy, LDR, inverse treatment planning, compressed sensing, real-time (Some figures may appear in colour only in the online journal) 1. Introduction Despite declining mortality rates (Center et al 2012), prostate carcinoma still accounts for more than 258 000 deaths worldwide and hence remains the second most common male cancer (Ferlay et al 2012). The majority of cases in the United States are diagnosed as localized disease (Brawley 2012). Local curative treatment options at this stage include radical prostatectomy (RPE), external beam radiation therapy (EBRT), and brachytherapy (Mohler et al 2010). Randomized comparisons in clinical trials between these treatment options have not been completed. The ongoing German PREFERE trial will address this issue but is not expected to be completed before 2030 (Stoeckle and Bussar-Maatz 2012). Nevertheless, comparative analyses indicate modern brachytherapy compares favorably in terms of both prostate-specific antigen (PSA) free survival and therapeutic side effects with RPE and EBRT (Grimm et al 2012). Brachytherapy treatment of prostate carcinoma was first described in 1909 (Aronowitz et al 2011). In 1981 Kumar and Bartone (1981) introduced the transperineal approach with iodine-125 seeds for continuous low dose rate (LDR) brachytherapy, which is preceding modern LDR brachytherapy with intra-operative dosimetric planning (Kaplan et al 2000). Over the last 20 years, different approaches for inverse brachytherapy treatment planning (ITP) have been proposed. The first approaches were semi-automatic (Bauer-Kirpes et al 1988, Roy et al 1991), while recent strategies apply stochastic global optimization strategies like simulated annealing (SA) (Pouliot et al 1996) or genetic algorithms (GA) (Yu 1996). SA and GA are general-purpose heuristics that are not tailored to the underlying global optimization problem structure. This allows using them for multiple applications without modification of the code. However, we often observe problem-specific heuristics outperforming these strategies. Very recently, the ITP was reformulated as a mixed-integer linear problem (MILP) (Yoo et al 2003, Yoo et al 2007). Using branch-and-bound with linear programming (LP) relaxation, which is one of the most efficient strategies for locating the global optimum for this type of problem structure (D'Souza et al 2001), we could show that global solutions can be located in reasonable time (Guthier and Hesser 2014). However, the solutions were found at a much higher cost than traditional SA and GA. The novelty of this paper is the reformulation of ITP in a mathematical structure that is of a compressed sensing (CS) type objective function. This allows using highly efficient matching pursuit like solvers that are generally considered among the fastest optimization strategies for CS problems (Du et al 2012). The research question addresses whether such performance gains can be observed for ITP as well, which in turn enables for the first time real-time planning updates while performing the intervention. The paper is structured as follows. The theory section focuses on the ITP and the underlying objective function. A derivation of the CS formulation of ITP follows. In material and methods, a matching pursuit inspired algorithm for inverse planning (MPIP) is introduced as a solver. Parameters and settings for the experimental study are explained. The result section describes the comparison of MPIP against the most efficient state-of-the-art solver so far on one hand side. On the other side, in order to estimate a lower bound of the objective function, a relaxed version of the ITP objective function is used to identify how near the solution is to this bound and hence the global optimum. Finally, MPIP is compared against a commercial TPS. Discussion and conclusion finish the paper. 2180

C Guthier et al

Phys. Med. Biol. 60 (2015) 2179

2. Theory 2.1.  LDR dose calculation and optimization

The LDR treatment planning goal is both to cover a planning target volume (PTV) with a userdefined lower dose bound and to spare organs at risks (OARs) by formulating organ-specific upper dose bounds. One approach of a dose-based objective function (DBOF) Q (x) penalizes values above and/or below given dose bounds (Baltas 2006). Defining a set of l potential seed positions S = {ξ1, ⋅, ξl }, ξi ∈ ℝ3, binary vector x ∈ {0, 1}l labels the occupancy of a position by a seed (value 1) or its absence (value 0). In the following, the objective function for the inverse planning process, Q (x), penalizes over and underdosages (relative to the given bounds) at given dose points for different organs. Specifically, if ν ∈ O = {1, ⋅, m} is the index set labeling considered organs and P ν = {ρν1, … , ρνn} is a set of three dimensional dose points ρνi ∈ ℝ3 inside the organ ν, we formulate Q  (x) =



ν∈O

1 ν R (x ) ; ∣P ν∣

(1)

Hereby, ∣Pν∣ is the total number of dose points of organ ν. The penalty term Rν (x) is the objective function related to organ ν and is defined as a weighted sum over violations of upper and lower bounds at each dose point: Rν (x) = w Lν 



ρ∈Pν

f Lν (x, ρ) + w Hν



ρ∈Pν

f Hν (x, ρ)

(2)

f Lν (x, ρ) penalizes a potential underdosage using the weighting factor w Lν. f Hν (x, ρ) accounts for overdosage with weighting factor w Hν . These penalty terms f Lν , f Hν are linear functions of the deviation from the given bound:

f Lν (x, ρ) = Θ {D Lν (ρ) − d(x, ρ) } · (D Lν (ρ) − d(x, ρ) )

(3)

f Hν (x, ρ) = Θ {d(x, ρ) − D Hν (ρ) } · (d(x, ρ) − D Hν (ρ) )

(4)

where Θ{η} is the Heaviside step function, d(x, ρ) is the dose received at the dose point ρ from a given seed configuration x. In addition, D Lν (ρ) is the lower and D Hν (ρ) the upper dose bound of the organ ν. For treatment planning, upper (D Hν (ρ) ) and lower (D Lν (ρ) ) dose bounds are assigned for PTV. Both address the clinical objective that PTV should be covered by the prescribed dose which is at least D Lν (ρ) and does not exceed D Hν (ρ). For OARs, their upper dose limit D Hν (ρ) is related to the maximum acceptable dose, whereas their lower dose limit is zero (and has hence not to be introduced in the objective function since the dose values are always positive). The dose d(x, ρ) is given by: ⎛ d(  x, ρ) = ⎜ ⎝

∫0

ttreat

. ⎞ d (r (ξ, ρ) , θ (ξ, ρ) , t ) ⎟ ·dt ·x ⎠

(5)

= (6)  : A·x . with d (r (ξ, ρ) , θ (ξ, ρ) , t ) is the dose-rate, r = r (ξ, ρ) is the Euclidean distance between seed ξ and dose point ρ, θ (ξ, ρ) is the angle between the vector (ξ − ρ) and the source long axis (Rivard et al 2004). In addition, ttreat is the overall treatment time. In the case of permanent LDR, it is set to infinity. For the dose-rate calculation, recommendations from AAPM TG-43 (Anderson 1995, Rivard et al 2004, Rivard et al 2007) are used. A, which is considered as a matrix Aij = d(r(ξi, ρj), 2181

C Guthier et al

Phys. Med. Biol. 60 (2015) 2179

θ(ξi, ρj)), is defined as dose dictionary where for each potential seed position a dose distribution is pre-calculated, an operation that is performed within ⩽ 0.2 seconds on a CPU. S containing all potential seed positions is constructed from a set of needles N  =  {1, ⋅, o}. In the discussed case of prostate brachytherapy, each needle trajectory is given by the transperineal template. This trajectory intersects with the PTV such that the potential seed positions are within an interval bounded by (χis, χie). In a general case, we have to consider multiple intervals per needle; which can be treated by the same formalism. Over this 1D interval [χis, χie] potential seed positions are evenly distributed, i.e.: ⎧ ⎧ p − 1⎫ ⎫ 1 ⎬⎬ ,…, i Γi = ⎨ξ∣ξ = λ χis + (1 − λ ) χie : λ = ⎨0, pi − 1 ⎭ ⎭ ⎩ pi − 1 ⎩ ⎪







(7)

where pi is the total number of potential seed positions on the respective needle. The set of all available seed positions S = i ∈ N Γi. Those dose points ρ are removed that are within a user-defined distance (in our case 1.5 mm) to the seeds i where xi = 1 (selected seeds). This distance is called the seed margin. To address the minimization of utilized needles as well a λ weighted penalty term on the number of used needles N (x) is used, reading: ∼  (x)= Q(x)+ λ · N(x). (8) Q



DBOF is often referred to as a dose volume histogram- (DVH-) based objective function where the dose bounds formulate the DVH that is to be achieved in the optimal case (Baltas 2006). 2.2.  Inverse treatment planning as compressive sensing problem

The goal of ITP is to minimize DBOF: x * = argmin min I x ∈ {0,1}

Q(x)

(9)

In appendix A, we show that this optimization problem can be reformulated as a MILP problem. As supported by general experience in brachytherapy ITP (Guthier and Hesser 2014), we assume in the remainder of this paper a sparse set of seeds in the global optimum for this problem. For reformulation, we ask for a plan with a guaranteed quality tolerance Q0 that minimizes the number of seeds used, i.e.: x * = argmin min I ‖ x ‖0 x ∈ {0,1}

s. t.

Q(x) ≤ Q0

(10)

where ∥x∥0 = {# i ∈ I ∣ xi ≠ 0 } is the l0 semi-norm, which counts the number of the nonzero elements of x. x* is hereby the sparse solution of the problem. The above formulated problem has the same mathematical structure as those for compressed sensing and is henceforth called in the following compressed sensing problem (CSP) (see appendix A). Over the last ten years, the performance of CSP solvers has improved by orders of magnitude. Among the best strategies are Matching Pursuit (Mallat and Zhang 1993) (MP), Orthogonal Matching Pursuit (Pati et al 1993) (OMP), Compressive Sampling Matching Pursuit (Needell and Tropp 2008) (CoSaMP), Subspace Pursuit (Dai and Milenkovic 2009) (SP), and variations thereof (Elad 2009). The underlying strategy sequentially chooses elements from the dictionary A to construct the solution. Although such greedy algorithms do not guarantee identifying the global optimum, under certain requirements on the dictionary elements the global optimum can be found with high probability (Tropp 2004). Whether our problem fulfills these requirements is beyond the scope of this paper. 2182

C Guthier et al

Phys. Med. Biol. 60 (2015) 2179

3.  Materials and methods 3.1.  Algorithm for CS in ITP

ITP formulated as a compressed sensing problem enables us to use tailored heuristics for optimization of the objective function. Matching pursuit (MP), as the basic strategy of a class of efficient compressed sensing heuristics, constructs a solution by iteratively including those seeds, elements of the dictionary A, which reduce the objective function Q (.) most. The iteration stops when the defined quality tolerance Q0 is reached or the objective function value cannot be decreased further. As main drawback of this method, the number of seeds grows continuously, leading to suboptimal solutions with respect to sparsity and DBOF. SP or CoSaMP overcomes this limitation. Both methods are able to identify unsuited components and exclude them from the solution found so far. In comparison, SP as well as CoSAMP outperforms MP in sparsity and quality of the solution (Needell and Tropp 2008, Dai and Milenkovic 2009). To fit to our problem we implemented the above strategies as follows: The algorithm starts with an empty support S (0) = 0 (i.e. no seeds) and the initial solution vector is hence x(0) = 0. During each iteration k the algorithm performs the following expansion and reduction step. ALGORITHM 1: MPIP strategy for solving the problem: vector x(k) represents the

solution at iteration k where each entry xi(k ) represents a potential seed position of the set S. The support S (k ) is the set containing indices of non-zero entries in x(k) and Q (x(k)) is the DBOF value. At each iteration step the algorithm selects the seed that minimizes the objective function value most and it’s index is added to the support S (k ). Afterwards, it is checked to see if a better solution can be achieved by removing one of the seeds of S (k ) and if so, x and S (k ) are updated correspondingly. If the DBOF value is below the given bound Q0 or the objective function cannot be improved anymore the iteration is stopped. Algorithm: MPIP Task: Approximate the solution of the ITP-problem (equation (9)) Definition: ⎧ ⎫ ⎧ 1 if i = j ei(j ) = ⎨       e(j ); , i, j =⎨1, … , l ⎬ i j   ≠ 0 if ⎩ ⎩ ⎭ Initialization:   x (o) = 0, S (o) = 0 Iteration:   do    k = k + 1

S ̂ ={1, ⋯ , m}

∼     j* = arg  minj ∈ S \̂ S (k − 1)Q (x (k − 1) + e(j ))     x(k ) = x(k − 1) + e(j* ),

S (k ) = S (k − 1) ∪ {j* } ∼ (k )     j* = arg min j ∈ S (k − 1)Q (x − e(j )) ∼ ∼     if   Q (x(k ) − e(j* ))

A new optimization method using a compressed sensing inspired solver for real-time LDR-brachytherapy treatment planning.

This work discusses a novel strategy for inverse planning in low dose rate brachytherapy. It applies the idea of compressed sensing to the problem of ...
964KB Sizes 1 Downloads 8 Views