Sparsity constrained split feasibility for dose-volume constraints in inverse planning of intensity-modulated photon or proton therapy
S. Penfold R. Zalas M. Casiraghi M. Brooke Y. Censor R. Schulte





Sparsity Constrained Split Feasibility for Dose-Volume Constraints in
Inverse Planning of Intensity-Modulated Photon or Proton Therapy
    Scott Penfold^1,2† The contributions of the first two authors to
this work are of equal shares., Rafał Zalas^3†, Margherita
Casiraghi^4
Mark Brooke^2, Yair Censor^5, Reinhard Schulte^6
^1Department of Medical Physics, Royal Adelaide Hospital
Adelaide, SA 5000, Australia 
 ^2Department of Physics, 
University of Adelaide, Adelaide, SA 5005, Australia
^3Department of Mathematics, The Technion, Technion City,
Haifa 32000, Israel 

^4Paul Scherrer Institute, Center for Proton Therapy (CPT)
Switzerland 

^5Department of Mathematics, University of Haifa
Mt. Carmel, Haifa 3498838, Israel 

^6Department of Basic Sciences, School of Medicine 
Loma Linda University, Loma Linda, CA 92354, USA.

(scott.penfold@sa.gov.au)

    July 10, 2016; Revised January 10, 2017
===================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




A split feasibility formulation for the inverse problem of intensity-modulated radiation therapy (IMRT) treatment planning with dose-volume constraints (DVCs) included in the planning algorithm is presented. It involves a new type of sparsity constraint that enables the inclusion of a percentage-violation constraint in the model problem and its handling by continuous (as opposed to integer) methods. We propose an iterative algorithmic framework for solving such a problem by applying the feasibility-seeking CQ-algorithm of Byrne combined with the automatic relaxation method (ARM) that uses cyclic projections. Detailed implementation instructions are furnished. Functionality of the algorithm was demonstrated through the creation of an intensity-modulated proton therapy plan for a simple 2D C-shaped geometry and also for a realistic base-of-skull chordoma treatment site. Monte Carlo simulations of proton pencil beams of varying energy were conducted to obtain dose distributions for the 2D test case. A research release of the Pinnacle^3 proton treatment planning system was used to extract pencil beam doses for a clinical base-of-skull chordoma case. In both cases the beamlet doses were calculated to satisfy dose-volume constraints according to our new algorithm. Examination of the dose-volume histograms following inverse planning with our algorithm demonstrated that it performed as intended. The application of our proposed algorithm to dose-volume constraint inverse planning was successfully demonstrated. Comparison with optimized dose distributions from the research release of the Pinnacle^3 treatment planning system showed the algorithm could achieve equivalent or superior results.


Keywords: dose-volume constraints, intensity-modulated radiation therapy, sparsity constraints, split feasibility, the CQ-algorithm, inverse planning, automatic relaxation method.





§ INTRODUCTION 


Intensity-modulated radiation therapy (IMRT) with photons or intensity-modulated proton therapy (IMPT) are rapidly evolving techniques for planning and delivering radiation therapy to solid tumors. For many tumor sites, IMRT with photons has superseded standard radiation therapy (RT) techniques and is becoming the new standard in RT delivery <cit.>. At existing proton centers, IMPT in combination with active pencil beam scanning is increasingly being used, replacing older passively scattered and collimated proton therapy techniques as a means for more accurately delivering high doses to the target volume and sparing of organs at risk (OARs) as indicated by dosimetric studies <cit.>. 

Instead of using a single upper dose bound for OARs and single lower dose bound for the target volumes, it has become a common practice in clinical trials and off-trial photon IMRT treatments to specify more than one dose-volume constraint (DVC), allowing a certain percentage of volume to violate to a certain extent a given bound. This additional DVC, which could be single or multiple, rely on accumulated clinical experience with conformal RT techniques. For example, Gulliford et al. <cit.> performed a detailed dose-volume analysis of the incidence of clinically relevant late rectal toxicities in patients treated with high-dose photon IMRT for prostate cancer and found that the incidence of moderate-to-severe rectal toxicity for any of six late-toxicity endpoints decreased incrementally for patients whose treatment plans met increasing numbers of DVCs from the set of V_30≤80%, V_40≤65%, V_50≤55%, V_60≤40%, V_65≤30%, V_70≤15%, and V_75≤3%. Here, V_X ≤ Y corresponds to a dose-volume constraint that Y% of the volume cannot receive more than X Gy. These and similar DVCs for OARs have found their way into clinical trial protocols and practice guidelines over the years, see, e.g., <cit.>.

Most modern inverse planning algorithms attempt to incorporate DVCs by defining sub-volumes with different dose objectives applied to each sub-volume. The multiple objectives are then combined into a single cost function to be minimized. Minimization in RT inverse planning with DVCs has been performed with a number of different approaches. Spirou and Chui <cit.> used gradient descent to seek a vector of ray intensities that minimized a cost function representing the sum of all dose constraints violations. However, incorporating DVCs directly into the cost function of the minimization process often renders the objective function non-convex and non-differentiable. This has the disadvantage of potentially resulting in local minima and thereby sub-optimal treatment plans. Cho et al. <cit.> used a similar concept but applied simulated annealing for minimization. Simulated annealing is less susceptible to non-convexity and non-differentiability but is less computationally efficient than gradient descent. Romeijn et al. <cit.> adopted a linear programming approach to handle what they called partial-volume constraints. However to make the problem tractable for computation, they replaced the familiar concept of DVCs by a closely related, but not identical notion of conditional value-at-risk (C-VaR). Zhang and Merritt <cit.> proposed a new least-squares model to handle DVCs while retaining differentiability at the expense of having to deal with a nested double minimization problem. Therefore, an inverse planning algorithm for DVCs that is computationally efficient, robust to non-convexity and non-differentiability yet without simplifying the problem statement has yet to be developed.



In the current work, feasibility-seeking methods, as opposed to minimization algorithms, are applied to RT inverse planning with DVCs. Within the proposed feasibility-seeking approach issues of convexity and differentiability of the cost function do not arise at all because no cost function is used. While the DVCs do require a constraint that is not convex (the sparsity-norm constraint set), we are able to incorporate it into the projection method that we use to solve the feasibility-seeking problem. This is possible because we have devised a way to calculate the projection onto this set in spite of it being non-convex.





Another general advantage of the feasibility-seeking approach has to do with the availability of a class of highly efficacious feasibility-seeking projection methods. These methods refer to iterative algorithms that use projections onto sets while relying on the general principle that when a family of, usually closed and convex, constraints sets is present, then projections onto the individual sets are easier to perform than projections onto other sets (intersections, image sets under some transformation, etc.) that are derived from the individual sets. Furthermore, projection methods may have algorithmic structures that are particularly suited for parallel computing, such as block-iterative projections (BIP) or string-averaging projections (SAP). They also demonstrate desirable convergence properties and good initial behavior patterns. See, for example, the 1996 review <cit.>, the recent annotated bibliography of books and reviews <cit.> and its references, and <cit.>.









We recently showed that IMPT inverse planning is possible with a fully-discretized, feasibility-seeking approach by iteratively projecting solution vectors in the beam intensity vector space onto half-spaces representing dose constraints in target and OAR volumes <cit.>. In our preliminary work, we demonstrated that with these iterative projection algorithms, feasible solutions meeting the planning objectives can be found that meet target and normal tissues dose bounds, in particular, if the constraints are not too challenging and/or the treatment modality is very conformal (e.g., by using protons).

In this paper, we use the fully-discretized feasibility-seeking approach applicable to either photon IMRT or IMPT inverse planning which leads to a mathematical feasibility problem. The upper and lower bounds on the doses to the various structures define the linear inequality constraints of the feasibility problem, which is solved by feasibility-seeking projection methods without attempting to minimize any cost function. Within this setup, we propose and investigate a novel method for allowing the feasibility-seeking inverse planning algorithm to automatically account for DVCs.

In the next section, we rigorously define the notion of percentage-violation constraint (PVC), which does not seem to have been used in the mathematical optimization community until now. A PVC injects integers into the problem which makes it difficult to solve. To circumvent this difficulty, we reformulate the PVC with the aid of a sparsity norm that counts the number of non-zero entries in a vector. This enables us to replace the original feasibility problem with PVC by another feasibility problem that includes non-convex constraints for the sparsity norm. For the resulting feasibility problem with this non-convex sparsity norm induced constraint we develop a new iterative projection algorithm which is a combination of the CQ-algorithm <cit.> and the automatic relaxation method (ARM) <cit.>.




§ METHODS



 §.§ Linear feasibility with percentage-violation constraints


Given p closed convex subsets Q_1,Q_2,⋯,Q_p⊆ R^n of
the n-dimensional Euclidean space R^n, expressed as level sets

    Q_j={  x∈ R^n| f_j(x)≤ v_j}  , for all
    j∈ J:={1,2,…,p},

where f_j:R^n→ R are convex functions and v_j are some
given real numbers, the convex feasibility problem (CFP) is to find a point
x^∗∈ Q:=∩_j∈ JQ_j. If Q=∅ where ∅ is the empty set then the CFP is said to
be inconsistent.

Denoting the inner product of two vectors in R^n by ⟨ a,b⟩:=∑_i=1^na_i b_i, we consider the following linear feasibility problem (LFP) with
percentage-violation constraint (PVC).


Linear Feasibility with Percentage-Violation
Constraint (PVC). Given a CFP as in (<ref>) with f_j(x)= ⟨
a^j,x⟩ and two real numbers 0≤α≤1 and
0<β<1, find a vector x^∗ that solves the system

    ⟨ a^j,x⟩≤ v_j, for all j∈ J

subject to the additional PVC constraint that:

    𝐈𝐧 𝐮𝐩 𝐭𝐨 𝐚𝐟𝐫𝐚𝐜𝐭𝐢𝐨𝐧 α(𝐢.𝐞.,100α%)𝐨𝐟 𝐭𝐡𝐞 𝐭𝐨𝐭𝐚𝐥
    𝐧𝐮𝐦𝐛𝐞𝐫𝐨𝐟𝐢𝐧𝐞𝐪𝐮𝐚𝐥𝐢𝐭𝐢𝐞𝐬
    
       𝐢𝐧 (<ref>) 𝐭𝐡𝐞 𝐫𝐢𝐠𝐡𝐭
    -𝐡𝐚𝐧𝐝 𝐬𝐢𝐝𝐞 𝐛𝐨𝐮𝐧𝐝𝐬v_j
     𝐦𝐚𝐲 𝐛𝐞 𝐩𝐨𝐭𝐞𝐧𝐭𝐢𝐚𝐥𝐥𝐲 𝐯𝐢𝐨𝐥𝐚𝐭𝐞𝐝
       𝐛𝐲 𝐮𝐩 𝐭𝐨 𝐚𝐟𝐫𝐚𝐜𝐭𝐢𝐨𝐧 β
    (𝐢.𝐞.,100β%)𝐨𝐟𝐭𝐡𝐞𝐢𝐫𝐯𝐚𝐥𝐮𝐞𝐬.




A PVC is an integer constraint by its nature. It changes the CFP to which it
is attached from being a continuous feasibility problem into becoming a mixed
integer feasibility problem. In the field of intensity-modulated radiation therapy (IMRT) treatment
planning dose-volume constraints (DVCs) are traditionally used to evaluate
treatment plans. DVCs are percentage-violation constraints but without
properly incorporating them into the algorithm itself it is not possible to a
priori guarantee that a solution will indeed obey them.

In this paper we propose a novel way to incorporate PVCs via the notion of a sparsity
norm and derive a tractable model and algorithmic approach, along with
detailed implementation instructions for using it, to solve DVCs feasibility
problems for inverse planning in IMRT.



 §.§ IMRT problem statement 


We consider the following linear interval feasibility problem (LIFP) which is
the basic model for the inverse problem in the fully-discretized approach to
IMRT treatment planning <cit.>:


Linear Interval Feasibility: the basic model
for the inverse problem in the fully-discretized approach to IMRT treatment
planning. Find x^∗∈ R^n for which the following hold:

    0≤ A_1x≤ b^1,
    
    b^3≥ A_2x≥ b^2,
    
    0≤ A_3x≤ b^4,
    
    x≥0,

where A_1∈ R_+^m_1× n, A_2∈ R_+^m_2× n,
A_3∈ R_+^m_3× n are given matrices, b^1∈ R_+^m_1,
b^2,b^3∈ R_+^m_2, b^4∈ R_+^m_3 are given vectors.
(The subscript + denotes the nonnegative orthant.)


In IMRT the row inequalities of (<ref>) represent voxels of an organ at
risk (OAR) whose permitted absorbed doses should not exceed b_t^1 for
each voxel t in this structure. The row inequalities of (<ref>) represent voxels of another OAR whose
permitted absorbed doses should not exceed b_t^4 for each voxel t in
this structure. The row inequalities of (<ref>) represent voxels of a planning
target volume (PTV) whose permitted absorbed doses should be above
b_t^2, but should not exceed b_t^3, for each voxel t in this
structure.

Our tool to “translate” the integer
constraint (<ref>) into a “continuous” one is the notion of sparsity norm, called
elsewhere the zero-norm, of a vector x∈ R^n which counts the number of
nonzero entries of x, that is,

    ‖ x‖_0:=|{x_i| x_i≠0}| ,

where |·| denotes here the cardinality, i.e., the number of elements
of a set. This notion has been recently used for various purposes in
compressed sensing, machine learning and more. The “lower +
 operation” on a vector x∈ R^n means that, for all
i=1,2,…,n,

    (x_+)_i:=max(0,x_i)={[ x_i, if x_i>0,;   0, if x_i≤0. ].  .

Obviously, x_+ is always a component-wise nonnegative vector. Hence,
‖ x_+‖_0 counts the number of positive entries of x and is
defined by

    ‖ x_+‖_0:=|{x_i| x_i>0}|.


To incorporate a DVC related to (<ref>) into the LIFP of Problem
<ref> we formulate another feasibility problem as follows.


Linear Interval Feasibility with DVC for the
inverse problem in the fully-discretized approach to IMRT treatment planning.
Find x^∗∈ R^n for which

    0≤ A_1x≤(1+β)b^1,
    
    b^3≥ A_2x≥ b^2,
    
    0≤ A_3x≤ b^4,
    
    x≥0,
    ‖(A_1x-b^1)_+‖_0≤α m_1,

where A_1, A_2, A_3, b^1, b^2, b^3 and b^4 are as
in (<ref>)–(<ref>), and β>0 and α∈0,1] are
given real numbers.


In this problem (<ref>) allows the doses to voxels of this
structure to “overflow” by β.
(<ref>) represents an OAR to which we do not attach a DVC for now.
(<ref>) represents a PTV to which we do not attach a DVC for now.
(<ref>) are the nonnegativity constraints on the solution vector of intensities.

The novelty of the model lies in (<ref>). It says that since we demanded
originally A_1x≤ b^1 in (<ref>) we must look at the
“plussed difference vector” (A_1
x-b^1)_+. It is nonnegative and has a nonzero component exactly and only
in components that belong to row inequalities in (<ref>) for which
(<ref>) is violated.

The zero-norm of (A_1x-b^1)_+ is thus equal to the number of those
violations and (<ref>) restricts this number to be not greater than
α m_1 where m_1 is the total number of row inequalities (i.e.,
voxels) in the OAR described by (<ref>). Thus, (<ref>)
guarantees that the number of violations up to β in (<ref>)
remains at bay under the number α m_1. In the following we
propose to use an efficient iterative projections method to solve Problem
<ref>.



 §.§ Projection methods for feasibility-seeking


Projections onto sets are used in a wide variety of methods in optimization
theory but here projection methods refer to iterative algorithms that use
projections onto sets while relying on the general principle that when a
family of, usually closed and convex, sets is present, then projections onto
the given individual sets are easier to perform than projections onto other
sets (intersections, image sets under some transformation, etc.) that are
derived from the given family of individual sets.

Projection methods may have different algorithmic structures, such as
block-iterative projections (BIP) or string-averaging projections (SAP) of
which some are particularly suitable for parallel computing, and they
demonstrate nice convergence properties and/or good initial convergence
patterns. This class of algorithms has witnessed great progress in recent
years and its member algorithms have been applied with success to many
scientific, technological and mathematical problems. See, e.g., the 1996
review <cit.>, the recent annotated bibliography of books and reviews
<cit.> and its references, the excellent book <cit.>, or
<cit.>.

For the LIFP of Problem <ref> one can use any of a variety of
projection methods to handle linear inequality constraints. The most famous of
those might be the Agmon-Motzkin-Schoenberg (AMS) cyclic feasibility-seeking
algorithm <cit.>. In this paper we adopt a projection
method of a particular nature, namely, the automatic relaxation method (ARM)
for solving interval linear inequalities of <cit.>.

ARM has two advantages over other projection methods applicable to this
problem: (i) it handles in each iteration an interval constraint and does not
need to handle the right-hand side and left-hand side inequalities of an
interval separately, (ii) additionally, it automatically implements a
relaxation strategy for the projections which takes into account how far from
the hyperslab, defined by an interval constraint, is the point that needs to
be projected on it and automatically and continuously adjusts the relaxation
parameter for the projection accordingly. The ARM generalizes the algebraic
reconstruction technique ART3 <cit.> and is further discussed
in Subsection 5.10 of Censor and Zenios <cit.>.



 §.§ Algorithmic approach


First we observe that Problem <ref> is a split feasibility
problem. Split feasibility
problems were introduced first in <cit.> and further studied in
<cit.> and many other publications. The constraints (<ref>)–(<ref>) can be collectively
described by c≤ Ax≤ b, where A is an (m_1+m_2+m_3)× n
matrix composed from blocks

    A:=(
    c]l
    A_1
    A_2
    A_3
     )  ,

b is an (m_1+m_2+m_3) vector given by

    b:=(
    c]l
    (1+β)b^1
    b^3
    b^4
     )  ,

and c is an (m_1+m_2+m_3) vector given by

    c:=(
    c]l
    0
    b^2
    0
     )  ,

and they, along with (<ref>) all reside in the space R^n of
intensity vectors x. On the other hand, the sparsity constraint
(<ref>) takes place in the space R^m_1 where the vectors of
doses in the OAR (<ref>) are, namely, the vector b^1 and the vectors
y=A_1x. Therefore, we must use not plain feasibility-seeking methods but
feasibility-seeking methods for split feasibility problems.

In the space R^n of intensity vectors we define the set

    C:={x∈ R^n| c≤ Ax≤ b}∩ R_+^n

where A, b and c are as in (<ref>), (<ref>) and (<ref>),
respectively, and R_+^n is the nonnegative orthant of R^n. In
R^m_1, the space of dose vectors of the OAR structure represented by
(<ref>), we define the set

    Q:={y∈ R^m_1|‖(y-b^1)_+‖_0≤α m_1}

with b^1 and α m_1 as in (<ref>). If a point y=A_1x
is in Q then it is guaranteed to fulfil (<ref>). So, our split
feasibility problem is to find a point x^∗∈ C such that A_1
x^∗∈ Q, precisely describing Problem <ref> above.

Common feasibility or split feasibility problems deal with convex sets
but here we observe that Q is not a convex set. However, we show below how
to project onto it orthogonally, thus enabling to use a
feasibility-seeking projection method for our Problem <ref>.

To solve the split feasibility formulation of Problem <ref> we propose to use the CQ-algorithm <cit.> for the sets C and Q
given by (<ref>) and (<ref>), respectively. It has the advantage
that it does not require to calculate the inverse A_1^-1 of A_1 in
order to “go back” from R^m_1 to
R^n within the iterative process. Instead, it uses the transposed
matrix A_1^T which is readily available. The CQ-algorithm
<cit.> is in fact a projected Landweber method for
the split feasibility formulation of Problem <ref>.

In the sequel P_Ω(z) denotes an orthogonal projection of a vector z
onto a set Ω. All data quantities mentioned below are as in Problem
<ref>. Since Q is not a convex set there might be more than one point for P_Q in (<ref>) below, therefore, the symbol ∈ therein means that x^k+1 could be any projection point onto Q of the vector in the parentheses whose projection onto Q is sought after, and can be
arbitrarily chosen from those if more then one exists.



Next we explain how to do the projections onto C and onto Q, and how to choose the parameter γ in (<ref>). Since Q of (<ref>) is not convex, the projection P_Q may by
multivalued. Nevertheless, for any z∈ R^m_1, we can calculate
P_Q(z) by using the following formula


    P_Q(z)=P_Q(z-b^1)+b^1

where

    Q:={y∈ R^m_1|‖ y_+‖_0≤α m_1}.

Hence the projection of a point z∈ R^m_1 onto the set Q of
(<ref>) is obtained by projecting the shifted point (z-b^1) onto the
set Q and adding b^1 to the result. The proof of this fact can be
found in the Appendix.

Therefore, the problem reduces to computing a projection onto Q. This is done according to the following recipe: First count how many
components of (z-b^1)
are positive, say ℓ. Then,

    P_Q(z-b^1)={[    (z-b^1), If ℓ≤α m_1,;          w, If ℓ>α m_1, ].

where w is the vector obtained from (z-b^1) by replacing its
ℓ-α m_1 smallest positive components by zeros and leaving the
others unchanged. If ℓ≤α m_1 then the point (z-b^1) is
already inside Q, thus P_Q(z-b^1)=(z-b^1). We will use
the above for z= A_1x^k in
(<ref>).

Following the seminal CQ-algorithm <cit.>, designed for the
case when both sets C and Q are convex, we propose that the parameter
γ in (<ref>) will be user-chosen from the open
interval 0<γ<2/θ where θ is pre-calculated once. To do so we employ <cit.> by using the
squared Frobenius matrix norm ‖ A_1‖_F^2 and defining

    θ:=‖ A_1‖_F^2=∑_i=1^m_1
    
    
    ∑_j=1^n
    | a_ij| ^2,

where for i=1,2,…,m_1 and j=1,2,…,n, the entries of A_1
are a_ij.

In the practical implementation we replace the projection onto C
(<ref>) by a sequence of projections onto the individual
inequalities of the constraints (<ref>)–(<ref>) that are
collectively described by c≤ Ax≤ b with where A, b and c are as
in (<ref>), (<ref>) and (<ref>), respectively, according to a
feasibility-seeking projection method of our choice. All of the above leads to our proposed Dose-Volume Split-feasibility (DVSF) Algorithm.



Algorithm <ref> is a general scheme that is made specific by choosing a feasibility-seeking projection
method to be used in its Step 3. Consult Bauschke and Borwein <cit.> for a review of such algorithms, see Censor and Cegielski <cit.> for an annotated bibliography of books and reviews on the subject and Censor et al. <cit.> for a review with experimental results.

We adopted here the automatic relaxation method (ARM) for feasibility-seeking <cit.>. We give
a generic description of this algorithm by considering the problem of solving
iteratively large and possibly sparse systems of interval linear inequalities
of the form

    w_j ≤ ⟨ a^j,x⟩≤ v_j,   j=1,2,...,p,

where a^j∈ R^n are given, for all j, and w=(w_j)∈ R^p, and
v=(v_j)∈ R^p are given too. Assuming that the system is
feasible, an x^∗∈ R^n which solves (<ref>) is
required. Geometrically, the system represents p nonempty hyperslabs in
R^n, each being the nonempty intersection of a pair of half-spaces. If we
are willing to ignore the slabs structure of the problem it could be addressed
as a system of 2p linear one-sided inequalities and solved by the Agmon-Motzkin-Schoenberg (AMS)
algorithm <cit.>. The ARM takes advantage of the interval structure of the problem by
handling in every iterative step a pair of inequalities and it also
realizes a specific relaxation principle (see <cit.> for details) in an
automatic manner. External relaxation parameters are available on top of the
built-in automatic relaxation principle.

For every hyperslab of the system (<ref>) denote by

    H_j:={x∈ R^n | ⟨ a^j,x⟩
    =v_j} and H_j:={x∈ R^n | ⟨
    a^j,x⟩=w_j}

its bounding hyperplanes. The median hyperplane will be

    H_j:={x∈ R^n | ⟨ a^j,x⟩=1/2(v_j+w_j)},

and the half-width ψ_j of the hyperslab is

    ψ_j=v_j-w_j/2∥ a^j∥,

where ∥ . ∥ stands for the Euclidean 2-norm.
The signed distance of a point z∈ R^n from the j-th median hyperplane
H_j is given by

    d(z,H_j)=⟨ a^j,z⟩-1/2(v_j+w_j
    )/∥ a^j∥.



Denoting d_j(k):=d(x^k,H_j(k)), the automatic relaxation method is as follows.





 §.§ Performance Testing

Performance tests with two different geometries were carried out to verify the functionality of the proposed algorithmic structure for IMRT. Applications to IMPT are presented in the current work. However, the algorithm is not proton specific, and is equally applicable to any form of IMRT. Only the values of the matrix A differ when different forms of radiation are used.



  §.§.§ Simplified 2D C-shaped geometry

A 2D test geometry was defined to simulate an axial cross-section of a tumour volume surrounding an organ at risk. The test geometry is illustrated in Figure <ref>. Structure pixels were defined with a resolution of 1 mm, also coinciding with the dose grid.



A proton pencil beamlet spacing of 2 mm, evenly distributed throughout the PTV structure, was used. Three beam angles were used to deliver dose to the PTV area. Each beam contained 146 proton pencil beamlets. The dose deposited by each pencil beamlet in the dose grid was calculated with the Monte Carlo toolkit Geant4 <cit.> and recorded in a text file.

The simulated beamlets were uniform circular proton beams of 2 mm
diameter. A pre-absorber made of 5.5 cm of polyethylene was inserted in the
beams at 5 cm in front of the irradiated geometry in order to smooth the Bragg
peaks and avoid dose distribution ripples due to beamlet spacing. The
beamlet energies for each aiming point were extracted from a calibration curve.
The energy used ranged from 118.5 MeV to 153 MeV with a resolution of 0.5 MeV.
The material of all the structures of the irradiated geometry was assumed to be
water.

The standard electromagnetic physics (G4EmStandardPhysics) and hadron physics
models (G4HadronPhysicsQGSP_BIC_HP) were used for proton tracking. Hadron
elastic scattering physics, stopping physics, ion physics and decays models were
also activated. A range cut of 0.1 mm was set for all particles. For each beamlet, 10^6 events were simulated and the mean absorbed dose per proton was calculated at each pixel of the dose grid.

A series of dose-volume constraints (DVCs) were defined to verify the functionality of the algorithm. These included:


	
  * dose only constraints (DOCs) applied to both the PTV and OAR structures
	
  * a single DVC associated with a single structure (the OAR structure)
	
  * multiple (two) DVCs associated with a single structure (the OAR structure)
	
  * DVCs associated with multiple structures (the PTV structure and the OAR structure)


At this point it is instructive to reconcile the dose-volume terminology used in the current work and the terminology commonly used in the literature. Let us consider an example where a prescription has been made to an OAR such that only 20% of the volume can receive more than 40 Gy and none of the volume can receive more than 50 Gy. Using the terminology of the current work, this would correspond to values of α = 0.2, b_1 = 40 Gy, and β = 0.25 in Problem <ref>. Using the common terminology, this would correspond to D_20%≤ 40 Gy and D_max = 50 Gy.

Table <ref> lists the combination of DVCs enforced in the current work, using the common terminology of IMRT DVCs. The dose prescriptions and percentage volume violations were chosen to allow for a demonstration of the functionality of the algorithm.

The initial pencil beam intensity vector before inverse planning was set to unity. The dose distribution resulting from the initial intensity vector is shown in Figure <ref>(a). The proposed algorithm was run for 2000 cycles for each prescription listed in Table <ref>. In this terminology, one cycle corresponds to one complete processing of all DVCs and DOCs applied to each pixel within both the PTV and OAR structures.





  §.§.§ Clinical 3D geometry

In keeping with the 2D geometry, a base of skull chordoma IMPT treatment plan was chosen due to the challenging constraints imposed by a target structure surrounding an avoidance structure. The Philips Pinnacle^3 treatment planning system (Philips Healthcare, Koninklijke Philips N.V.) was used to contour the PTV and brainstem. The exported DICOM RT (structure) files were imported into a MATLAB (The MathWorks, Inc.) script and the brainstem and PTV contours were mapped over the CT coordinates. A dose grid was created in MATLAB to match that defined in Pinnacle^3. The dimensions were 42 × 43 × 9 voxels with resolutions of 2 mm, 2 mm and 3 mm in the x, y and z dimensions, respectively. The dose grid was twice as large as the CT pixel size in the x and y dimension and equivalent to the CT resolution in the z dimension. A reduced number of slices (9) was required due to memory restrictions encountered during the export of pencil beamlet doses.

An IMPT treatment plan was created in the Pinnacle^3 research release of proton pencil beam scanning (PBS). Two beams were targeted at the PTV from angles of 80^o and 280^o, containing 574 and 564 beamlets, respectively. A range shifter of 7.5 cm thickness was used with both beams to ensure proximal PTV coverage. Distal and proximal margins for pencil beam placement were automatically calculated as a percentage of proton range. The dose grid resulting from each unit intensity beamlet was exported from Pinnacle^3. Beamlet parameters were set to 80% layer overlap, a lateral spot resolution of 0.6 cm, a lateral target margin of 0.4 cm and 3 standard deviation dose spread during dose calculation. Dose was calculated with the analytical PBS algorithm which includes nuclear attenuation and an energy and material dependent multiple Coulomb scattering model.



For each structure A-matrices were created by combining the geometry defined by the DICOM RT structures and the dose grid obtained for each beamlet. Each 3D beamlet dose grid was rearranged to a 1D vector which became a column of an A-matrix. Each row of the A_OAR matrix corresponded to a voxel of the brainstem and likewise each row of the A_PTV matrix corresponded to a voxel of the PTV.

Two DVCs were tested for the base of skull chordoma IMPT treatment plan (see Table <ref>). The DVCs differed in the dose objectives for the brainstem while keeping the PTV objectives constant. The same DVCs were applied consistently for both the DVSF algorithm (Algorithm <ref>) and Pinnacle^3.



Independent values for the parameter γ of (<ref>) were used for the OAR and PTV and are denoted by γ_PTV and γ_OAR. These values were determined from the structure-specific calculation of θ in (<ref>), denoted by θ_PTV and θ_OAR. The relaxation parameters λ_k of (<ref>) are fixed throughout the iterations and represented by λ_PTV and λ_OAR.



§ RESULTS



 §.§ Simplified 2D C-shaped geometry

The dose distributions following inverse planning for Prescriptions 1 and 4 in Table <ref> are shown in Figure <ref>(b) and <ref>(c). The dose-volume histograms following inverse planning for all cases listed in Table <ref> are presented in Figure <ref>.



The dose distributions (Figure <ref>) allow for a qualitative assessment of the functionality of the DVSF algorithm (Algorithm <ref>). It is evident that the dose resulting from unit intensity pencil beamlets is successfully modulated toward the desired dose distribution. However, for a quantitative assessment the dose-volume histograms must be considered. When Prescription 1 DOCs were applied the dose objectives on the PTV structure could not be met (Figure <ref>(a)). Introducing the DVC on the OAR structure relaxed these conditions and resulted in satisfaction of the dose objectives on the PTV structure (Figure <ref>(b)). While the DVC on the OAR structure was not achieved in Prescription 2, continued iterations would have resulted in a dose distribution approaching the DVC more closely. The DVSF algorithm (Algorithm <ref>) was shown to function with multiple DVCs applied to a single structure (Figure <ref>(c)), and with DVCs applied to multiple structures (Figure <ref>(d)).





 §.§ Clinical 3D geometry

Cumulative DVHs for Prescription 1 of Table <ref> using the DVSF algorithm (Algorithm <ref>) and that produced by the Pinnacle^3 inverse planning algorithm are shown in Figure <ref>. All constraints of the less challenging dose objectives were met by the DVSF algorithm (Algorithm <ref>) whereas Pinnacle^3 exceeded the maximum dose for the OAR and did not satisfy the PTV minimum dose constraint. It should be noted that the Pinnacle^3 inverse planning was run only once with unit weighting on all dose objectives. It is possible that alteration of the objective weightings by trial-and-error may have resulted in a more desirable dose distribution. However, the objective of the current work was to compare the inherent ability of the algorithms to satisfy the inverse problem, and as such, iterative plan refinement by altering objective weights was not considered. Dose distributions for Prescription 1 of Table <ref> are shown in a single axial slice in Figure <ref>. The DVSF algorithm (Algorithm <ref>) showed higher conformality of the target structure.



Cumulative DVHs for Prescription 2 of Table <ref> are shown in Figure <ref>(b). It is clear that both the DVSF algorithm (Algorithm <ref>) and Pinnacle^3 had more difficulty meeting the dose objectives in this case. The DVSF algorithm (Algorithm <ref>) was better able to meet the hard dose constraints when compared to Pinnacle^3 but the latter was closer to meeting the D_95≥ 70 Gy DVC on the PTV. Dose distributions for Prescription 2 of Table <ref> are shown in a single axial slice in Figure <ref>. Both dose distributions show cold spots in the target region.







§ DISCUSSION AND CONCLUSION

A new DVSF algorithm (Algorithm <ref>) based on feasibility-seeking has been successfully applied to IMPT inverse planning in the current work. The proposed DVSF algorithm (Algorithm <ref>) is based on a modification of the CQ-algorithm of Byrne <cit.> and is capable of directly incorporating the DVCs associated with radiation therapy prescriptions into the split feasibility-seeking problem statement. Our DVSF algorithm (Algorithm <ref>) is not restricted to IMPT and is equally applicable to other forms of IMRT inverse planning. Test cases consisted of a simplified 2D C-shape target surrounding an avoidance structure and a clinical base of skull chordoma abutting the brainstem.

The DVSF algorithm (Algorithm <ref>) performs orthogonal projections to satisfy both the DVCs and the lower and upper dose constraints. The AMS cyclic projection method <cit.> was implemented for single-sided inequality dose objectives and the ARM algorithm of <cit.> was implemented for interval inequalities (i.e., upper and lower dose bounds for a given structure).

A series of experiments were performed with 2D C-shaped geometry using varying DVCs to validate the functionality of our DVSF algorithm (Algorithm <ref>). While DVC aims were not met in all cases within the allowed number of iterations, the shape of the DVH curve verified that the algorithm was attempting to meet these objectives. Experimentation with user-defined relaxation parameter values γ and λ was performed to investigate the effect of these settings on algorithmic performance. When λ was left at the fixed value of 1, it was found that γ values closer to the upper allowable limit of 2/θ were required to meet the DVC aims. Further work concerning automatic choice of these user-defined parameters is currently being undertaken and will be presented in future investigations.

A clinical 3D IMPT treatment geometry was also investigated. The performance of the DVSF algorithm (Algorithm <ref>) was compared to that of the research release of Pinnacle^3 with proton pencil beam scanning. The shape of DVHs differed for the two inverse planning algorithms. For the prescriptions investigated, our DVSF algorithm (Algorithm <ref>) was found to result in a more conformal dose distribution when assessing isodose contours and DVH distributions. It is acknowledged that the dose distributions obtained with Pinnacle^3 may be improved with the addition of planning structures. However, to allow for a comparison of the inverse planning algorithms directly, no such structures were included in the treatment planning method.

While the implementation of the DVSF algorithm (Algorithm <ref>) was sequential in the current work, the structure of the algorithm lends itself to parallelization. For example, block-iterative or string-averaging projection operators may be used when performing the orthogonal projections described in Step 3 of Algorithm <ref>. Such implementations will not only have benefits in computational speed, but may also result in superior dose distributions, as has been observed in the use of these algorithms in tomographic image reconstruction <cit.>. Further work will examine the potential of block-iterative and string-averaging algorithmic schemes for the DVSF algorithm (Algorithm <ref>).



§ APPENDIX


Here is a proof of formula (<ref>) for the projection calculation onto
the non-convex set Q.


We show that the following translation formula

    P_Q(z)=P_Q(z-b^1)+b^1

holds true for every z∈ R^m_1, despite the fact that P_Q and
P_Q are set-valued, i.e., a point z∈ R^m_1 might have more
than one projection onto the set. Note that

    Q=Q-b^1.

By the definition of projection of a point onto a set,

    q_0∈ P_Q(z) if and only if q_0∈ Q and ‖ z-q_0‖≤‖
    z-q‖, for all q∈ Q.

Similarly, (q_0-b^1)∈ P_Q(z-b^1) if and only if (q_0-b^1)∈Q̂ and

    ‖(z-b^1)-(q_0-b^1)‖≤‖(z-b^1)-q‖

holds for every q∈Q. Therefore, by (<ref>),
(<ref>) and (<ref>), we have the following
equivalences

    q_0∈ P_Q(z)     ⟺ q_0∈ Q and
     ∀_q∈ Q‖ z-q_0‖≤‖ z-q‖
       ⟺(q_0-b^1)∈Q and
       ∀_q∈ Q ‖(z-b^1)-(q_0-b^1)‖≤‖(z-b^1
    )-(q-b^1)‖
       ⟺(q_0-b^1)∈Q and
       ∀_q∈Q ‖(z-b^1)-(q_0-b^1)‖≤‖(z-b^1)-q‖
       ⟺(q_0-b^1)∈ P_Q(z-b^1)
       ⟺ q_0∈ P_Q(z-b^1)+b^1,

which completes the proof.


Acknowledgments. We thank the two anonymous referees for their constructive comments which helped us improve the paper. The work of Y. Censor and R. Schulte was supported by Research Grant No. 2013003 of the United States-Israel Binational Science Foundation (BSF) and by Award No. 1P20183640-01A1 of the National Cancer Institute (NCI) of the National Institutes of Health (NIH). The authors thank Philips for technical assistance with Pinnacle^3 software for this research. 
























































































































ieeetr


