Sparsity constrained split feasibility for dose-volume constraints in inverse planning of intensity-modulated photon or proton therapy

TL;DR

This paper introduces a novel sparsity-constrained split feasibility approach for inverse planning in intensity-modulated radiation therapy, effectively incorporating dose-volume constraints using continuous methods, demonstrated through proton therapy case studies.

Contribution

It presents a new split feasibility formulation with sparsity constraints for dose-volume constraints in IMRT planning, enabling continuous optimization methods and improved treatment plans.

Findings

Successfully applied to 2D and clinical proton therapy cases

Achieved dose distributions satisfying dose-volume constraints

Performed as well as or better than existing treatment planning systems

Abstract

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 Pinnacle3 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 Pinnacle3 treatment planning system showed the algorithm could achieve equivalent or superior results.

Tables

Table 1. Table 1: Prescriptions associated with PTV and OAR structures in order to test the functionality of the proposed DVSF algorithm (Algorithm 5 ) in a simplified 2D geometry.

Prescription	OAR	PTV
1	$D_{max}=45$ Gy	$D_{min}=70$ Gy $D_{max}=77$ Gy
2	$D_{15.5\%}\le 45$ Gy $D_{max}=70$ Gy	$D_{min}=70$ Gy $D_{max}=77$ Gy
3	$D_{50\%}\le 25$ Gy $D_{15.5\%}\le 45$ Gy $D_{max}=70$ Gy	$D_{min}=70$ Gy $D_{max}=77$ Gy
4	$D_{8.5\%}\le 45$ Gy $D_{max}=70$ Gy	$D_{min}=66.5$ Gy $D_{95\%}\ge 70$ Gy $D_{max}=77$ Gy

Table 2. Table 2: Prescriptions associated with PTV and OAR (brainstem) structure for a clinical test case.

Prescription	OAR	PTV
1	$D_{\text{max}}=54$ Gy	$D_{\text{min}}=66.5$ Gy
	$D_{5\%}\le 50$ Gy	$D_{\text{max}}=74.9$ Gy
		$D_{95\%}\ge 70$ Gy
2	$D_{\text{max}}=40$ Gy	$D_{\text{min}}=66.5$ Gy
	$D_{5\%}\le 35$ Gy	$D_{\text{max}}=74.9$ Gy
		$D_{95\%}\ge 70$ Gy

Equations

$$Q_j=\left\{x\in R^n\mid f_j(x)\le v_j\right\},\text{ for all }j\in J:=\{1,2,\dots ,p\},$$

$$⟨a^j,x⟩\le v_j,\text{ for all }j\in J$$

$$\mathbf{In}\mathbf{up}\mathbf{to}\mathbf{a}\mathbf{fraction}\alpha (i.e.,100\alpha \%)\mathbf{of}\mathbf{the}\mathbf{total}\mathbf{number}\mathbf{of}\mathbf{inequalities}$$

$$\mathbf{in}\text{ (refeq:original-ineq) }\mathbf{the}\mathbf{right}\text{-}\mathbf{hand}\mathbf{side}\mathbf{bounds}{v}_j\mathbf{may}\mathbf{be}\mathbf{potentially}\mathbf{violated}$$

$$\mathbf{by}\mathbf{up}\mathbf{to}\mathbf{a}\mathbf{fraction}\beta (i.e.,100\beta \%)\mathbf{of}\mathbf{their}values.$$

$$0\le A_{1}x\le b^1,$$

$$b^3\ge A_{2}x\ge b^2,$$

$$0\le A_{3}x\le b^4,$$

$$x\ge 0,$$

$$\Vert x{\Vert }_0:=\left|\{x_i\mid x_i\ne 0\}\right|,$$

$$(x_{+})_{i}:=\max(0,x_{i})=\left\{\begin{array}[c]{l}x_{i},\text{ if }x_{i}>0,\\ 0,\text{ if }x_{i}\leq 0.\end{array}\right..$$

$$\Vert x_{+}{\Vert }_0:=|\{x_i\mid x_i>0\}|.$$

$$0\le A_{1}x\le (1+\beta )b^1,$$

$$b^3\ge A_{2}x\ge b^2,$$

$$0\le A_{3}x\le b^4,$$

$$x\ge 0,$$

$$\Vert (A_{1}x-b^1{)}_{+}{\Vert }_0\le \alpha m_1,$$

$$A:=\left(\begin{tabular}[c]{l}$A_{1}$\\ $A_{2}$\\ $A_{3}$\end{tabular}\ \right),$$

$$b:=\left(\begin{tabular}[c]{l}$(1+\beta)b^{1}$\\ $b^{3}$\\ $b^{4}$\end{tabular}\ \right),$$

$$c:=\left(\begin{tabular}[c]{l}$0$\\ $b^{2}$\\ $0$\end{tabular}\ \right),$$

$$C:=\{x\in R^n\mid c\le Ax\le b\}\cap R_{+}^n$$

$$Q:=\{y\in R^{m_1}\mid \Vert (y-b^1{)}_{+}{\Vert }_0\le \alpha m_1\}$$

$$x^{k+1}\in P_C\left(x^k+\gamma A_1^T\left(P_Q(A_1{x}^k)-A_1{x}^k\right)\right).$$

$$P_Q(z)=P_{\overline{Q}}(z-b^1)+b^1$$

$$\overline{Q}:=\{y\in R^{m_1}\mid \Vert y_{+}{\Vert }_0\le \alpha m_1\}.$$

$$P_{\overline{Q}}(z-b^{1})=\left\{\begin{array}[c]{ll}(z-b^{1}),&\text{If }\ell\leq\alpha m_{1},\\ w,&\text{If }\ell>\alpha m_{1},\end{array}\right.$$

$$\theta :=\Vert A_1{\Vert }_F^2=\sum _{i=1}^{m_1}\sum _{j=1}^n{\left|a_{ij}\right|}^2,$$

$$v^k:=P_{\overline{Q}}(A_1{x}^k-b^1).$$

$$P_Q(A_1{x}^k)=v^k+b^1.$$

$$u^k=x^k+\gamma A_1^T\left(P_Q(A_1{x}^k)-A_1{x}^k\right).$$

$$w_j\le ⟨a^j,x⟩\le v_j,j=1,2,...,p,$$

$${\overline{H}}_j:=\{x\in R^n|⟨a^j,x⟩=v_j\}\text{ and }{\underline{H}}_j:=\{x\in R^n|⟨a^j,x⟩=w_j\}$$

$$H_j:=\{x\in R^n|⟨a^j,x⟩=\frac{1}{2}(v_j+w_j)\},$$

$${\psi }_j=\frac{v_j-w_j}{2\parallel a^j\parallel },$$

$$d(z,H_j)=\frac{⟨a^j,z⟩-\frac{1}{2}(v_j+w_j)}{\parallel a^j\parallel }.$$

Full text

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

1 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 [1, 2, 3, 4, 5]. 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 [6, 7, 8, 9, 10].

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. [11] 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}\leq 80\%$, $V_{40}\leq 65\%$, $V_{50}\leq 55\%$, $V_{60}\leq 40\%$, $V_{65}\leq 30\%$, $V_{70}\leq 15\%$, and $V_{75}\leq 3\%$. Here, $V_{X}\leq 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., [12].

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 [13, Section F] 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. [14] 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. [15, 16] 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 [17] 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 [18], the recent annotated bibliography of books and reviews [19] and its references, and [20].

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 [21]. 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 [22] and the automatic relaxation method (ARM) [23].

2 Methods

2.1 Linear feasibility with percentage-violation constraints

Given $p$ closed convex subsets $Q_{1},Q_{2},\cdots,Q_{p}\subseteq R^{n}$ of the $n$-dimensional Euclidean space $R^{n}$, expressed as level sets

[TABLE]

where $f_{j}:R^{n}\rightarrow R$ are convex functions and $v_{j}$ are some given real numbers, the convex feasibility problem (CFP) is to find a point $x^{\ast}\in Q:=\cap_{j\in J}Q_{j}.$ If $Q=\emptyset$ where $\emptyset$ is the empty set then the CFP is said to be inconsistent.

Denoting the inner product of two vectors in $R^{n}$ by $\left\langle a,b\right\rangle:=\sum_{i=1}^{n}{a_{i}b_{i}}$, we consider the following linear feasibility problem (LFP) with percentage-violation constraint (PVC).

Problem 1

Linear Feasibility with Percentage-Violation Constraint (PVC). Given a CFP as in (1) with $f_{j}(x)=$ $\left\langle a^{j},x\right\rangle$ and two real numbers $0\leq\alpha\leq 1$ and $0<\beta<1,$ find a vector $x^{\ast}$ that solves the system

[TABLE]

subject to the additional PVC constraint that:

[TABLE]

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.

2.2 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 [24, 25, 26]:

Problem 2

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

[TABLE]

where $A_{1}\in R_{+}^{m_{1}\times n}$, $A_{2}\in R_{+}^{m_{2}\times n}$, $A_{3}\in R_{+}^{m_{3}\times n}$ are given matrices, $b^{1}\in R_{+}^{m_{1}}$, $b^{2},b^{3}\in R_{+}^{m_{2}}$, $b^{4}\in R_{+}^{m_{3}}$ are given vectors. (The subscript + denotes the nonnegative orthant.)

In IMRT the row inequalities of (4) 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 (6) 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 (5) 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 (3) into a “continuous” one is the notion of sparsity norm, called elsewhere the zero-norm, of a vector $x\in R^{n}$ which counts the number of nonzero entries of $x,$ that is,

[TABLE]

where $\mid\cdot\mid$ 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\in R^{n}$ means that, for all $i=1,2,\ldots,n,$

[TABLE]

Obviously, $x_{+}$ is always a component-wise nonnegative vector. Hence, $\|x_{+}\|_{0}$ counts the number of positive entries of $x$ and is defined by

[TABLE]

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

Problem 3

Linear Interval Feasibility with DVC for the inverse problem in the fully-discretized approach to IMRT treatment planning. Find $x^{\ast}\in R^{n}$ for which

[TABLE]

where $A_{1}$, $A_{2}$, $A_{3}$, $b^{1}$, $b^{2}$, $b^{3}$ and $b^{4}$ are as in (4)–(7), and $\beta>0$ and $\alpha\in[0,1]$ are given real numbers.

In this problem (11) allows the doses to voxels of this structure to “overflow” by $\beta.$ (13) represents an OAR to which we do not attach a DVC for now. (12) represents a PTV to which we do not attach a DVC for now. (14) are the nonnegativity constraints on the solution vector of intensities.

The novelty of the model lies in (15). It says that since we demanded originally $A_{1}x\leq b^{1}$ in (4) 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 (11) for which (4) is violated.

The zero-norm of $(A_{1}x-b^{1})_{+}$ is thus equal to the number of those violations and (15) restricts this number to be not greater than $\alpha m_{1}$ where $m_{1}$ is the total number of row inequalities (i.e., voxels) in the OAR described by (11). Thus, (15) guarantees that the number of violations up to $\beta$ in (11) remains at bay under the number $\alpha m_{1}.$ In the following we propose to use an efficient iterative projections method to solve Problem 3.

2.3 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 [18], the recent annotated bibliography of books and reviews [19] and its references, the excellent book [27], or [20].

For the LIFP of Problem 3 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 [28, 29]. In this paper we adopt a projection method of a particular nature, namely, the automatic relaxation method (ARM) for solving interval linear inequalities of [23, Algorithm 1].

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 [30] and is further discussed in Subsection 5.10 of Censor and Zenios [31].

2.4 Algorithmic approach

First we observe that Problem 3 is a split feasibility problem. Split feasibility problems were introduced first in [32] and further studied in [33, 34] and many other publications. The constraints (11)–(13) can be collectively described by $c\leq Ax\leq b$, where $A$ is an $(m_{1}+m_{2}+m_{3})\times n$ matrix composed from blocks

[TABLE]

$b$ is an $(m_{1}+m_{2}+m_{3})$ vector given by

[TABLE]

and $c$ is an $(m_{1}+m_{2}+m_{3})$ vector given by

[TABLE]

and they, along with (14) all reside in the space $R^{n}$ of intensity vectors $x.$ On the other hand, the sparsity constraint (15) takes place in the space $R^{m_{1}}$ where the vectors of doses in the OAR (4) are, namely, the vector $b^{1}$ and the vectors $y=A_{1}x.$ 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

[TABLE]

where $A$, $b$ and $c$ are as in (16), (17) and (18), 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 (4), we define the set

[TABLE]

with $b^{1}$ and $\alpha m_{1}$ as in (11). If a point $y=A_{1}x$ is in $Q$ then it is guaranteed to fulfil (15). So, our split feasibility problem is to find a point $x^{\ast}\in C$ such that $A_{1}x^{\ast}\in Q,$ precisely describing Problem 3 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 3.

To solve the split feasibility formulation of Problem 3 we propose to use the $CQ$-algorithm [22] for the sets $C$ and $Q$ given by (19) and (20), 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 [22, Algorithm 1.1] is in fact a projected Landweber method for the split feasibility formulation of Problem 3.

In the sequel $P_{\Omega}(z)$ denotes an orthogonal projection of a vector $z$ onto a set $\Omega.$ All data quantities mentioned below are as in Problem 3. Since $Q$ is not a convex set there might be more than one point for $P_{Q}$ in (21) below, therefore, the symbol $\in$ 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.

Algorithm 4

**The $CQ$-Algorithm for the Split Linear Feasibility Problem with a DVC. **

Step 0: Take an arbitrary $x^{0}\in R^{n}$, and set $k=0$.

Step 1: For a current iterate $x^{k}$ calculate $A_{1}x^{k}$ and compute the next iterate by

[TABLE]

If a stopping criterion applies then stop, otherwise go back to the Step 1 with $k\Longleftarrow k+1$.

Next we explain how to do the projections onto $C$ and onto $Q,$ and how to choose the parameter $\gamma$ in (21). Since $Q$ of (20) is not convex, the projection $P_{Q}$ may by multivalued. Nevertheless, for any $z\in R^{m_{1}},$ we can calculate $P_{Q}(z)$ by using the following formula

[TABLE]

where

[TABLE]

Hence the projection of a point $z\in R^{m_{1}}$ onto the set $Q$ of (20) is obtained by projecting the shifted point $(z-b^{1})$ onto the set $\overline{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 $\overline{Q}$. This is done according to the following recipe: First count how many components of $(z-b^{1})$ are positive, say $\ell$. Then,

[TABLE]

where $w$ is the vector obtained from $(z-b^{1})$ by replacing its $\ell-\alpha m_{1}$ smallest positive components by zeros and leaving the others unchanged. If $\ell\leq\alpha m_{1}$ then the point $(z-b^{1})$ is already inside $\overline{Q},$ thus $P_{\overline{Q}}(z-b^{1})=(z-b^{1}).$ We will use the above for $z=$ $A_{1}x^{k}$ in (21).

Following the seminal $CQ$-algorithm [22], designed for the case when both sets $C$ and $Q$ are convex, we propose that the parameter $\gamma$ in (21) will be user-chosen from the open interval $0<\gamma<2/\theta$ where $\theta$ is pre-calculated once. To do so we employ [35, Corollary 2.3] by using the squared Frobenius matrix norm $\|A_{1}\|_{F}^{2}$ and defining

[TABLE]

where for $i=1,2,\ldots,m_{1}$ and $j=1,2,\ldots,n$, the entries of $A_{1}$ are $a_{ij}.$

In the practical implementation we replace the projection onto $C$ (21) by a sequence of projections onto the individual inequalities of the constraints (11)–(13) that are collectively described by $c\leq Ax\leq b$ with where $A$, $b$ and $c$ are as in (16), (17) and (18), 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 5

**The Dose-Volume Split-feasibility (DVSF) Algorithm.

**

Step (-1): Read all data from Problem 3 and calculate (once) the transposed matrix $A_{1}^{T},$ the value of $\theta$ according to (25), and choose a parameter $\gamma$ in the open interval $0<\gamma<2/\theta$.

Step 0: Take an arbitrary $x^{0}\in R_{+}^{n}$, and set $k:=0$.

Step 1: Project $A_{1}x^{k}$ onto $Q$ as follows:

Step 1.1: For the current iterate $x^{k}$ compute $A_{1}x^{k}$, count the coordinates of $(A_{1}x^{k}-b^{1})$ that are positive and denote their number by $\ell.$

Step 1.2: Calculate (using (24) with $z=$ $A_{1}x^{k}$)

[TABLE]

Step 1.3: Calculate a projection of $A_{1}x^{k}$ onto $Q$ (following (22)–(23)):

[TABLE]

Step 2: Calculate $u^{k}\in R^{n}$ by the formula

[TABLE]

Step 3: Instead of projecting $u^{k}$ onto $C$ as required in (21), use $u^{k}$ from Step 2 as an initial point and perform a sweep (or several sweeps) of a feasibility-seeking projection method for the inequalities of (11)–(14). When stopping this sweep (or several sweeps) take the resulting vector as the next iterate $x^{k+1}$.

Step 4: If a stopping criterion applies then stop, otherwise go back to the Step 1 with $k\Longleftarrow k+1$.

Algorithm 5 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 [18] for a review of such algorithms, see Censor and Cegielski [19] for an annotated bibliography of books and reviews on the subject and Censor et al. [20] for a review with experimental results.

We adopted here the automatic relaxation method (ARM) for feasibility-seeking [23]. 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

[TABLE]

where $a^{j}\in R^{n}$ are given, for all $j$, and $w=(w_{j})\in R^{p}$, and $v=(v_{j})\in R^{p}$ are given too. Assuming that the system is feasible, an $x^{\ast}\in R^{n}$ which solves (29) 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 [28, 29]. 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 [23] 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 (29) denote by

[TABLE]

its bounding hyperplanes. The median hyperplane will be

[TABLE]

and the half-width $\psi_{j}$ of the hyperslab is

[TABLE]

where $\parallel.\parallel$ stands for the Euclidean 2-norm. The signed distance of a point $z\in R^{n}$ from the $j$-th median hyperplane $H_{j}$ is given by

[TABLE]

Denoting $d_{j(k)}:=d(x^{k},H_{j(k)}),$ the automatic relaxation method is as follows.

Algorithm 6

The Automatic Relaxation Method (ARM).

Initialization: $x^{0}\in R^{n}$ is arbitrary.

Iterative step: Given a current iterate $x^{k}$ calculate the next iterate $x^{k+1}$ by

[TABLE]

Control: The control sequence $\{j(k)\}_{k=0}^{\infty}$ according to which hyperslabs are picked during iterations is cyclic on $\{1,2,....,m\}$, i.e., $j(k)=k\bmod m+1$.

Relaxation parameters: External relaxation parameters $\{\lambda_{k}\}_{k=0}^{\infty}$, available on top of the built-in automatic relaxation principle are confined, for all $k\;\geq\;0,$ to

[TABLE]

2.5 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.

2.5.1 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 1. 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 [36] 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 $\alpha=0.2$, $b_{1}$ = 40 Gy, and $\beta=0.25$ in Problem 3. Using the common terminology, this would correspond to $D_{20\%}\leq 40$ Gy and $D_{max}=50$ Gy.

Table 1 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 2(a). The proposed algorithm was run for 2000 cycles for each prescription listed in Table 1. 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.

2.5.2 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 Pinnacle3 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 Pinnacle3. The dimensions were $42\times 43\times 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 Pinnacle3 research release of proton pencil beam scanning (PBS). Two beams were targeted at the PTV from angles of 80o and 280o, 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 Pinnacle3. 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 2). 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 5) and Pinnacle3.

Independent values for the parameter $\gamma$ of (28) were used for the OAR and PTV and are denoted by $\gamma_{PTV}$ and $\gamma_{OAR}$. These values were determined from the structure-specific calculation of $\theta$ in (25), denoted by $\theta_{PTV}$ and $\theta_{OAR}$. The relaxation parameters $\lambda_{k}$ of (34) are fixed throughout the iterations and represented by $\lambda_{PTV}$ and $\lambda_{OAR}$.

3 Results

3.1 Simplified 2D C-shaped geometry

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

The dose distributions (Figure 2) allow for a qualitative assessment of the functionality of the DVSF algorithm (Algorithm 5). 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 3(a)). Introducing the DVC on the OAR structure relaxed these conditions and resulted in satisfaction of the dose objectives on the PTV structure (Figure 3(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 5) was shown to function with multiple DVCs applied to a single structure (Figure 3(c)), and with DVCs applied to multiple structures (Figure 3(d)).

3.2 Clinical 3D geometry

Cumulative DVHs for Prescription 1 of Table 2 using the DVSF algorithm (Algorithm 5) and that produced by the Pinnacle3 inverse planning algorithm are shown in Figure 4. All constraints of the less challenging dose objectives were met by the DVSF algorithm (Algorithm 5) whereas Pinnacle3 exceeded the maximum dose for the OAR and did not satisfy the PTV minimum dose constraint. It should be noted that the Pinnacle3 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 2 are shown in a single axial slice in Figure 5. The DVSF algorithm (Algorithm 5) showed higher conformality of the target structure.

Cumulative DVHs for Prescription 2 of Table 2 are shown in Figure 4(b). It is clear that both the DVSF algorithm (Algorithm 5) and Pinnacle3 had more difficulty meeting the dose objectives in this case. The DVSF algorithm (Algorithm 5) was better able to meet the hard dose constraints when compared to Pinnacle3 but the latter was closer to meeting the $D_{95}\geq 70$ Gy DVC on the PTV. Dose distributions for Prescription 2 of Table 2 are shown in a single axial slice in Figure 6. Both dose distributions show cold spots in the target region.

4 Discussion and Conclusion

A new DVSF algorithm (Algorithm 5) based on feasibility-seeking has been successfully applied to IMPT inverse planning in the current work. The proposed DVSF algorithm (Algorithm 5) is based on a modification of the $CQ$-algorithm of Byrne [22] and is capable of directly incorporating the DVCs associated with radiation therapy prescriptions into the split feasibility-seeking problem statement. Our DVSF algorithm (Algorithm 5) 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 5) performs orthogonal projections to satisfy both the DVCs and the lower and upper dose constraints. The AMS cyclic projection method [28] was implemented for single-sided inequality dose objectives and the ARM algorithm of [23] 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 5). 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 $\gamma$ and $\lambda$ was performed to investigate the effect of these settings on algorithmic performance. When $\lambda$ was left at the fixed value of 1, it was found that $\gamma$ values closer to the upper allowable limit of $2/\theta$ 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 5) was compared to that of the research release of Pinnacle3 with proton pencil beam scanning. The shape of DVHs differed for the two inverse planning algorithms. For the prescriptions investigated, our DVSF algorithm (Algorithm 5) 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 Pinnacle3 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 5) 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 5. 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 [37]. Further work will examine the potential of block-iterative and string-averaging algorithmic schemes for the DVSF algorithm (Algorithm 5).

5 Appendix

Here is a proof of formula (22) for the projection calculation onto the non-convex set $Q$.

Proof. We show that the following translation formula

[TABLE]

holds true for every $z\in R^{m_{1}}$, despite the fact that $P_{Q}$ and $P_{\overline{Q}}$ are set-valued, i.e., a point $z\in R^{m_{1}}$ might have more than one projection onto the set. Note that

[TABLE]

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

[TABLE]

Similarly, $(q_{0}-b^{1})\in P_{\overline{Q}}(z-b^{1})$ if and only if $(q_{0}-b^{1})\in\hat{Q}$ and

[TABLE]

holds for every $\overline{q}\in\overline{Q}$. Therefore, by (37), (38) and (39), we have the following equivalences

[TABLE]

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 Pinnacle3 software for this research.