# Intensity-Modulated Radiation Therapy Optimization for Acceptable and Remaining-One Unacceptable Dose-Volume and Mean-Dose Constraint Planning

Authors: Ryosei Nakada; Omar M. Abou Al-Ola; Tetsuya Yoshinaga
Journal: Computational and Mathematical Methods in Medicine
 (2020)

---

## Abstract

We give a novel approach for obtaining an intensity-modulated radiation therapy (IMRT) optimization solution based on the idea of continuous dynamical methods. The proposed method, which is an iterative algorithm derived from the discretization of a continuous-time dynamical system, can handle not only dose-volume but also mean-dose constraints directly in IMRT treatment planning. A theoretical proof for the convergence to an equilibrium corresponding to the desired IMRT planning is given by using the Lyapunov stability theorem. By introducing the concept of “acceptable,” which means the existence of a nonempty set of beam weights satisfying the given dose-volume and mean-dose constraints, and by using the proposed method for an acceptable IMRT planning, one can resolve the issue that the objective and evaluation are different in the conventional planning process. Moreover, in the case where the target planning is totally unacceptable and partly acceptable except for one group of dose constraints, we give a procedure that enables us to obtain a nearly optimal solution close to the desired solution for unacceptable planning. The performance of the proposed approach for an acceptable or unacceptable planning is confirmed through numerical experiments simulating a clinical setup.


---

## Body

## 1. Introduction

Intensity-modulated radiation therapy (IMRT) [1, 2], which is a type of external beam radiotherapy, is an advanced radiotherapy technique used to reduce the amount of normal tissue exposed to radiation in the treatment area. It also improves the ability to conform to the volume of treatment applied to concave tumor forms. The radiation delivery pattern is determined using computerized applications designed to perform optimization and therapeutic simulations. The radiation dose is adjusted by modifying the intensity of the beam in accordance with the shape of the tumor. The intensity of the radiation dose is raised near the tumor but is reduced or zero among adjacent normal tissues, i.e., a radiation field with a nonuniform (i.e., modulated) intensity is delivered. Compared with conventional nonmodulated radiotherapy techniques with beams of uniform intensity, IMRT leads to better tumor targeting and mitigation of side effects, which may increase the chance of curing the patient. Furthermore, recently, volumetric modulated arc therapy (VMAT) has been adopted in clinical practice. One of the features of VMAT is that it performs IMRT while rotating the specific range of gantry, which increases the flexibility of generating a theoretically highly conformal treatment plan [3–5].The calculation of nonuniform intensities based on the dose prescription in the planning target volume (PTV) and the neighboring critical or sensitive organs, called organs at risk (OARs), is called inverse planning. IMRT inverse treatment planning uses optimization techniques, with the objective function measuring the quality of a treatment plan, to form the dose distribution with the ability to generate concave dose distributions and provide a specific spacing for the sensitive normal structures. The dose-volume constraint (DVC) expressed as an inequality condition with a set of a volume percentage and prescription dose applied to a PTV or OAR is evaluated on specific, predefined subvolumes of the organ and restricts the relative volume of a structure that receives more or less than a particular threshold. The use of DVCs is a normal way to determine the objective and has been the standard way of evaluating the treatment in practice; however, they are generally handled in limited ways in current optimization algorithms. Namely, because the inequality with percentages is difficult to treat as an objective function in optimization, the conventional methods are mostly constructed by using gradient techniques such as Newton’s method and the conjugate gradient method of minimizing dose-based and biology-based objective functions [6–11], and therefore, one expects to obtain a feasible solution violating the DVCs in a minimal way. Due to the difference between the objective for optimization and the evaluation for the planning result, a trial-and-error approach was required in the DVC-based planning process [12], more specifically as follows: When the given DVC is too tight, planners cannot easily obtain the optimal solution and the goal of treatment planning is not achieved. On the other hand, when the given DVC is too loose, planners can obtain the optimal solution, but there is room for improvement in IMRT treatment planning. In order to obtain an achievable and high-quality IMRT treatment planning, planners are required to find optimal DVC settings for each clinical case. As a result, it takes a lot of time to perform trial and error. Furthermore, the quality of the treatment planning varies from planner to planner [13]; a useful approach is required that supports planners in obtaining high-quality treatment planning in a short time. For example, on the hardware side, there is an approach to perform faster calculations utilizing a graphics processing unit (GPU) [14–16]. On the software side, there is an approach that uses a knowledge-based treatment plan with low variability [17, 18] and uses automatic planning [19] and a multicriteria optimization framework [20, 21] to reduce the frequency of planner intervention.In this paper, we present an optimization method to handle not only dose-volume but also mean-dose constraintsdirectly in IMRT treatment planning, which is an iterative algorithm derived from the discretization of a continuous-time dynamical system [22–29] with a set of equilibrium points. A theoretical proof for the convergence to an equilibrium corresponding to the desired IMRT planning is given by using the Lyapunov stability theorem. By introducing the concept of “acceptable,” which means the existence of a nonempty set of beam weights satisfying the given dose-volume and mean-dose constraints, and by using the proposed method for an acceptable IMRT planning, one can resolve the issue that the objective and evaluation are different in the conventional planning process. The system of nonlinear differential equations defined in this paper has a different vector field from that of the previously proposed system [29]. The previous method for obtaining an acceptable solution has a disadvantage in that it requires a long calculation time because of the high numerical cost of integrating piecewise differential equations describing the hybrid dynamical system. To resolve the disadvantage, we conduct a numerical discretization with multiplicative calculus [30, 31]. Moreover, for achieving a high-precision IMRT plan, we extend the previous continuous-time system to allow mean-dose as well as dose-volume concepts.Now, we can get an acceptable solution for any totally acceptable IMRT treatment planning; however, we cannot obtain a nearly optimal solution close to the desired solution for every unacceptable planning. Therefore, as a strategy to search for such a solution, we present an achievement procedure. That is, when assuming that the target planning is totally unacceptable and partly acceptable except for one group of dose constraints, the objective is to get a solution located as close as possible to the remaining-one group unacceptable constraint while satisfying the partly acceptable situation. In considering the dynamics in the state spaceΩ, schematically shown in Figure 1, the first step of the procedure is to restrict the trajectory within the partly acceptable subspace Ae using the proposed optimizing system for the objective function Ue such that the corresponding inverse problem is partly acceptable by excluding a function Ur from the total objective function with an acceptable set Ar with Ae∩Ar=ϕ. Then, the next step is to move the state x to decrease Urx and Uex in early and later iterations, respectively, by using the trajectory of an optimizing system, which has a time-dependent regularization effect [32] making it possible to optimize multiple criteria consisting of Ue and Ur, with an initial value ae in Ae. The optimizing system can get a desired solution in Ae that is close to the optimal subspace Ar for the remaining-one constraint function Ur.Figure 1
Schematic diagram of strategy for searching for nearly optimal solution.This paper is organized as follows. In Section2, some preliminaries and definitions are provided. In Section 3, we propose a novel iterative method for obtaining an acceptable dose-volume and mean-dose constraint planning based on the continuous-time dynamical system, where the convergence to an acceptable set of solutions is theoretically guaranteed. A procedure is presented in Section 4 to get a nearly optimal solution for the remaining-one unacceptable constraint planning. In Sections 5 and 6, the performance of the proposed approach for an acceptable or unacceptable planning is confirmed through numerical experiments simulating a clinical setup. Finally, in Section 7, some concluding remarks are drawn.In the rest of this section, we introduce notations that will be used below:θk indicates the kth element of the vector θ; and Θi and Θj indicate the ith row and jth column vectors of the matrix Θ, respectively.


## 2. IMRT Planning

The radiation oncologist performs IMRT planning to determine the PTV to be irradiated and the OARs to be spared. With IMRT, letJ be the total number of beamlets that have split off from the radiation beam, and this is calculated as the number of delivered beams multiplied by the number of beamlets in each beam head. The problem of dose calculation, in matrix notation, has the form
(1)D=Kx,where D∈R+I with R+ being the set of nonnegative real numbers which is the dose vector whose components represent the total dose deposited in each voxel of the patient’s three-dimensional volume, each component of x∈Ω¯⊂R+J with Ω¯ indicating the closure of the open hypercube Ω=0,γJ for some γ>0 which represents the intensity in the corresponding beamlet, and K∈R+I×J which consists of all elements kij representing the dose deposited in the ith voxel due to the unit intensity in the jth beamlet; it is called the dose influence (or deposition) matrix. If there are volumes of PTV and OAR with I1 and I2 voxels, respectively, then D includes D1∈R+I1 and D2∈R+I2 as subvectors, and K includes K1∈R+I1×J and K2∈R+I2×J as submatrices. A similar definition can be applied to the existence of multiple PTVs and OARs.Letb1L and b2U represent the lower and upper bounds on the dose delivered to PTV and OAR, respectively (b1L>0 and b2U≥0). Additionally, we define an upper dose bound on the dose delivered to PTV as b1U>b1L to avoid excessively high doses inside the PTV. We also handle a mean dose denoted by Dm as a real-valued function of the dose Dm∈R+Im, which is defined by
(2)Dm≔1Im∑i=1ImDmi,for m=1,2. With respect to the mean-dose constraints, the lower and upper bounds on the mean dose of radiation delivered to PTV and OAR are indicated by c1L and c2U, respectively.Definition 1.
The IMRT planning systemK1,b1U,K1,b1L,K2,b2U,K1,c1L,K2,c2U is consistent if the set
(3)E=e∈Ω:b1U≥K1ei1≥b1L,∀i1∈1,2,⋯,I1,K2ei2≤b2U,∀i2∈1,2,⋯,I2,K1e≥c1L,K2e≤c2Uis not empty; otherwise, it is inconsistent.Definition 2.
For the givenx∈Ω and each set of dose volumes and their conditions K1,b1U,ζ1U, K1,b1L,ζ1L, and K2,b2U,ζ2U, where ζ1U≤1, ζ1L≤1, and ζ2U≤1 are the prescribed proportion rates, the corresponding dose distribution is partly dose-volume acceptable if each number of elements of the index sets
(4)I1Ux=i1∈1,2,⋯,I1:K1xi1≤b1U,I1Lx=i1∈1,2,⋯,I1:K1xi1≥b1L,I2Ux=i2∈1,2,⋯,I2:K2xi2≤b2Uis, respectively, greater than the prescribed proportion of I1, I1, and I2, namely, each of the inequalities
(5)I1Ux≥ζ1UI1,(6)I1Lx≥ζ1LI1,(7)I2Ux≥ζ2UI2is satisfied for some x, where · indicates the number of elements in the set.Definition 3.
For the givenx∈Ω and each set of doses and their conditions (K1, c1L) and (K2, c2U), the corresponding dose distribution is partly mean-dose acceptable if each of the inequalities
(8)K1x≥c1L,K2x≤c2Uis satisfied.
We define that the IMRT planning system is acceptable if there exists a common beam set such that the dose distributions in PTVs and OARs are partly dose-volume and mean-dose acceptable for all dose-volume and mean-dose constraints.Definition 4.
The IMRT planning systemK1,b1U,ζ1U,K1,b1L,ζ1L,K2,b2U,ζ2U,K1,c1L,K2,c2U is acceptable if the following set is not empty:
(9)A=a∈Ω:I1Ua≥ζ1UI1,I1La≥ζ1LI1,I2Ua≥ζ2UI2,K1a≥c1L,K2a≤c2U.
If the IMRT planning system is consistent, then it is acceptable. We are interested in the situation where the system is inconsistent and acceptable. In this paper, the problem of dose-volume and mean-dose constraint optimization in IMRT planning is defined to obtain the unknown variablex∈A if the system is acceptable.
For describing the system, let us define the functionsP1UD1, P1LD1, P2UD2, Q1LD1, and Q2UD2 as follows:
(10)P1UD1i1=D1i1,ifD1i1≤b1U,b1U,otherwise,P1LD1i1=D1i1,ifD1i1≥b1L,b1L,otherwise,P2UD2i2=D2i2,ifD2i2≤b2U,b2U,otherwise,Q1LD1i1=D1i1,ifD1≥c1L,c1LD1D1i1,otherwise,Q2UD2i2=D2i2,ifD2≤c2U,c2UD2D2i2,otherwise,for i1=1,2,⋯,I1 and i2=1,2,⋯,I2. We consider the case where
(11)D≔D1D1D2D1D2,K≔K1K1K2K1K2,with I=3I1+2I2. To simplify the description, we also define
(12)PD≔P1UD1P1LD1P2UD2Q1LD1Q2UD2.


## 3. IMRT Optimization for Acceptable Dose-Volume and Mean-Dose Constraint Planning

This section provides the definition of the dose-volume and mean-dose constraint optimization method and its theoretical results. Consider an initial-value problem for the nonlinear differential equation(13)dxjtdt=xjt1−xjtγlogfjxtσj,where
(14)fjx=∑i1=1I1K1ji1P1UK1xi1K1xi1R1Ux+∑i1=1I1K1ji1P1LK1xi1K1xi1R1Lx+∑i2=1I2K2ji2P2UK2xi2K2xi2R2Ux+∑i1=1I1K1ji1Q1LK1xi1K1xi1S1Lx+∑i2=1I2K2ji2Q2UK2xi2K2xi2S2Ux,σj=3∑i1=1I1K1ji1+2∑i2=1I2K2ji2,for j=1,2,⋯,J with x0=x0∈Ω. Where
(15)R1Ux=0,ifxwithK1,b1U,ζ1Uispartlydose‐volumeacceptable,1,otherwise,R1Lx=0,ifxwithK1,b1L,ζ1Lispartlydose‐volumeacceptable,1,otherwise,R2Ux=0,ifxwithK2,b2U,ζ2Uispartlydose‐volumeacceptable,1,otherwise,S1Lx=0,ifxwithK1,c1Lispartlymean‐doseacceptable,1,otherwise,S2Ux=0,ifxwithK2,c2Uispartlymean‐doseacceptable,1,otherwise,for i1=1,2,⋯,I1 and i2=1,2,⋯,I2.We give theoretical results for the behavior of the solution to the dynamical system in Equation (13). First, we show that all solutions stay inside the hypercube.Proposition 5.
If we choose the initial valuex0∈Ω of the dynamical system in Equation (13), then the solution φt,x0 stays in Ω for all t>0.Proof.
Since the system can be written asdxj/dt=xj1−γ−1xjhjx with a function hj, we see that for any j, the solution satisfies dφj/dt≡0 on the subspace where xj=0 or xj=γ. Therefore, the subspace is invariant, and trajectories cannot pass through every invariant subspace in accordance with the uniqueness of solutions for the initial-value problem. This leads to any solution φt,x0 of the system in Equation (13) with initial value x0∈Ω being in Ω for all t>0.Next, we prove the stability of equilibrium in the setA, which corresponds to the desired radiation beam weights. Namely, the existence of a Lyapunov function for the system in Equation (13) guarantees the stability of the equilibrium set [33].Theorem 6.
If the IMRT planning system is acceptable, thenA in Equation (9) as an equilibrium set of Equation (13) is stable.Proof.
Any pointa∈A is an equilibrium of Equation (13), so we can say that A is an equilibrium set of equilibrium points. Consider a Lyapunov-candidate function defined in the set Ω as
(16)Vx=∑j=1Jσjajlogajxj+γ−ajlogγ−ajγ−xj=∑j=1Jσj∫ajxjv−ajv+v−ajγ−vdv=∑j=1Jσj∫ajxjγv−ajvγ−vdv,which is positive definite with respect to the point a∈A. We obtain its derivative along the solution to the system in Equation (13). When x∈Ω is neither partly dose-volume nor mean-dose acceptable, the derivative is given by
(17)V̇x13=∑j=1Jσjγxj−ajxjγ−xjẋj=−∑j=1Jajσjlogσj−1∑i=1IKjiPKxiKxi+∑j=1Jxjσjlogσj−1∑i=1IKjiPKxiKxi≤−∑j=1Jajσjσj−1∑i=1IKjilogPKxi−logKxi+∑j=1Jxjσjσj−1∑i=1IKjiPKxiKxi−1≤−∑i=1IKailogPKxi−logKxi+∑i=1IKxiPKxiKxi−1≤−∑i=1IKai1−KxiPKxi+∑i=1IKxiPKxiKxi−1=−∑i=1I1PKxiKxi−PKxiPKxi−Kai.
To treat anyx∈Ω, define the index sets:
(18)I¯1Ux=ϕ,ifxwithK1,b1U,ζ1Uispartlydose‐volumeacceptable,1,2,⋯,I1\I1Ux,otherwise,I¯1Lx=ϕ,ifxwithK1,b1L,ζ1Lispartlydose‐volumeacceptable,1,2,⋯,I1\I1Lx,otherwise,I¯2Ux=ϕ,ifxwithK2,b2U,ζ2Uispartlydose‐volumeacceptable,1,2,⋯,I2\I2Ux,otherwise,O¯1Lx=ϕ,ifxwithK1,c1Lispartlymean‐doseacceptable,1,2,⋯,I1,otherwise,O¯2Ux=ϕ,ifxwithK2,c2Uispartlymean‐doseacceptable,1,2,⋯,I2,otherwise.
Then, the calculation of the derivative is reduced to the following:(19)V̇x13=−1b1U∑i1∈I¯1UxK1xi1−b1Ub1U−K1ai1−1b1L∑i1∈I¯1Lxb1L−K1xi1K1ai1−b1L−1b2U∑i2∈I¯2UxK2xi2−b2Ub2U−K2ai2−1c1Lc1L−K1x∑i1∈O¯1LxK1ai1−c1LK1xK1xi1−1c2UK2x−c2U∑i2∈O¯2Uxc2UK2xK2xi2−K2ai2≤0.
The last inequality is supported by(20)∑i1=1I1K1ai1−c1LK1xK1xi1=I1K1a−c1L>0,∑i2=1I2c2UK2xK2xi2−K2ai2=I2c2U−K2a>0,for partly mean-dose unacceptable (K1, c1L) and (K2, c2U). The derivative is zero at x=a∈Ω. Thus, Vx is a Lyapunov function with respect to a. Consequently, the equilibrium set A is stable.A numerical discretization of differential equations describing the system in Equation (13) derives an iterative method. By using geometric multiplicative Euler discretization [30, 31], we obtain the iterative algorithm of the variable z∈Ω for radiation beam weight
(21)zjn+1=zjnexp1−zjnγlogfjznσjh,for j=1,2,⋯,J and n=0,1,2,⋯ with z0=x0∈Ω, where h>0 denotes a step size. When γ⟶∞, h=1, and all binary functions R1U, R1L, R1U, S1, and S2U are assumed to be one for simplicity, the jth element of z can be described as
(22)zjn+1=zjn∑i=1IKji∑i=1IKjiPKzniKzni,for j=1,2,⋯,J. We see that the iterative formula in Equation (22) is similar to that of the expectation-maximization (EM) method in computed tomography. Namely, replacing PKzn with the measured projection and assuming K and z are the projection operator and pixel density of the image to be reconstructed, respectively, the iterative formula describes the EM-type reconstruction algorithm.One can design a different combination of vector field and Lyapunov function from that of Equations (13) and (16). For example, a continuous analog of the iterative gradient algorithm based on the split feasibility problem [11] for handling dose constraints is able to be a continuous-time system with an equilibrium whose stability is supported by the Lyapunov theorem, but it is required to choose an appropriate value of the step size in the iterative formula as an additive Euler discretization for avoiding slow convergence and oscillatory behavior in solutions. Our resultant iterative algorithm as a modification of the EM procedure, which is a popular iterative image reconstruction method in practice, is expected to calculate the desired solution with a good convergence rate [34], when h=1.


## 4. IMRT Optimization for Remaining-One Unacceptable Dose-Volume and Mean-Dose Constraint Planning

Consider an IMRT planning system(23)S≔K1,b1U,ζ1U,K1,b1L,ζ1L,K2,b2U,ζ2U,K1,c1L,K2,c2U,K3,b3U,ζ3U,K3,c3U,where the subsystem
(24)K3,b3U,ζ3U,K3,c3U≕Sr,indicates an additional set of dose-volume and mean-dose constraints with the upper bound for an OAR, as an example of remaining-one unacceptable constraint, by rearranging the formulation described in Section 3. We assume that the planning system is totally unacceptable and the subsystem
(25)Se≔S\Sr=K1,b1U,ζ1U,K1,b1L,ζ1L,K2,b2U,ζ2U,K1,c1L,K2,c2U,is acceptable. We show a procedure to get a solution located as close as possible to the remaining-one unacceptable constraint of Sr while satisfying the partly acceptable situation as in the subsystem Se. In the procedure, we use an initial value problem for the following difference equation of variables zn with initial state zN1∈Ω starting at N1≥0:
(26)zjn+1=zjnexp1−zjnγlogfjzn+gjznσj+δjh,for n=N1,N1+1,⋯, where the functions gj and δj are, respectively, defined by
(27)gjz=∑i3=1I3K3ji3P3UK3zi3K3zi3R3Uz+∑i3=1I3K3ji3Q3UK3zi3K3zi3S3Uz,δj=2∑i3=1I3K3ji3,for j=1,2,⋯,J. Where
(28)P3UD3i3=D3i3,ifD3i3≤b3U,b3U,otherwise,Q3UD3i3=D3i3,ifD3≤c3U,c3UD3D3i3,otherwise,for i3=1,2,⋯,I3 and
(29)R3Uzn=0,ifzwithK3,b3U,ζ3Uispartlydose‐volumeacceptable,λn−N1,otherwise,S3Uzn=0,ifzwithK3,c3Uispartlymean‐doseacceptable,λn−N1,otherwise,with a positive parameter λ. Note that one takes λn⟶0 and then has gjzn⟶δj, j=1,2,⋯,J, for n⟶∞, when λ<1. We give the procedure as follows:Procedure 7.
For an IMRT planning systemS consisting of Equation (23), which is totally unacceptable and has a partly acceptable subsystem Se in Equation (25), the following shows a procedure for obtaining a solution located as close as possible to the remaining-one unacceptable constraint of Sr in Equation (24) while satisfying the partly acceptable situation as in the subsystem Se.Step 1.
The first step of the procedure is to solve iterative solutions to the difference equation in Equation (21) for the subsystem Se and obtain a stable fixed point ae in the set Ae illustrated in Figure 1. Obtaining ae∈Ae is guaranteed by Theorem 6 for the continuous analog of the difference system in Equation (21).Step 2.
The next step is to examine the iterative algorithm in Equation (26) with some value of the parameter λ<1 and an initial value zN1=ae∈Ae and move the state zN1+n, n=1,2,⋯,N2, while expecting the decrease of the objective function UrzN1+n, as shown in Figure 1. Note that because the partly acceptable situation of the subsystem Se means that fjz=σj in Equation (26) at z=ae∈Ae for j=1,2,⋯,J, the iterative steps can forcibly decrease the sequence Urznn=N1N1+ν for at least a small value of ν. By introducing the function gjzn with a time-varying parameter λn−N1 in Equation (26), the iterative process has a time-dependent regularization effect making it possible to optimize multiple criteria. As the function tends to δj independently of zn along with time n that has passed, the optimizing system can get a desired solution in Ae that is close to the optimal subspace Ar with respect to the objective function Ur for the remaining-one constraint.


## 5. Experimental Method

To evaluate the proposed method, we examine two IMRT treatment planning problems: “acceptable planning” and “remaining-one unacceptable planning.” Figure2 shows 128×128-sized phantom images for the planning study. The image in Figure 2(a) shows a phantom consisting of a core (blue region) and C-shape PTV (red region), which refers to a C-shape phantom published in the task group 119 report [35] from the American Association of Physicists in Medicine (AAPM). The image in Figure 2(b) represents a phantom simulating a head and neck with a PTV, parotid gland, and spinal cord, which are colored red, green, and blue, respectively. The former image was used for an acceptable planning, and the latter was used for a partly unacceptable planning. In many clinical cases, an ideal solution to an IMRT planning problem is not achieved; therefore, planning for IMRT requires a trial-and-error process to find a solution for “partly acceptable planning.” Generally, it is necessary for a planner to set priorities for each of the PTVs and OARs. For example, in the case of a simple structure phantom, such as the C-shape phantom, an ideal and feasible solution is obtained (“acceptable planning”) because the PTV and core can be handled with equal priority. However, in the case of a complex structure phantom, such as the head and neck phantoms, there is a case where the DVCs of other organs cannot be satisfied by achieving the DVCs of a high-priority organ. In this case, a planner prioritizes the PTV and spinal cord, and the parotid’s DVC is set to the next. In other words, the IMRT planning that satisfies the DVCs of the PTV and cord perfectly achieves the parotid’s DVC as much as possible. This is valuable because if the spinal cord is irradiated excessively, the risk of myelopathy caused by radiation will increase. In addition, to reduce the dose applied to the parotid gland as much as possible, the function of the salivary glands can be preserved [36–42].Figure 2
Phantom images of (a) C-shape for acceptable planning and (b) extended C-shape simulating the head and neck for remaining-one unacceptable planning. The red region shows PTV, and the blue and green regions indicate OARs.
(a)(b)Table1 shows the DVCs for the C-shape phantom case, and we adopted a 9-field beam arrangement with every 40-degree beam from 0 to 360 degrees, which is given as a harder version of the problem in the task group 119 report [35]. Table 2 shows the DVCs and mean doses for the head and neck phantom case. To achieve more stringent DVCs, we adopted a 36-field beam arrangement with every 5-degree beam from 0 to 175 degrees. Although there is research related to beam angle optimization in IMRT [43], the objective of the beam arrangement is to obtain dose distributions equivalent or superior to static gantry IMRT [3–5]. So we adopted a VMAT technique to perform complex treatment planning.Table 1
The prescribed constraints and equivalent parameters for acceptable planning (Vd denotes the percentage of volume receiving at least d [Gy]).
Assigned region (colored region in Figure2(a))OrganConstraintEquivalent parameterPTV (red)PTVV55<10%(K1,b1U,ζ1U) with b1U=55 and ζ1U=0.90V50≥95%(K1,b1L,ζ1L) with b1L=50 and ζ1L=0.95OAR (blue)CoreV10<5%(K2,b2U,ζ2U) with b2U=10 and ζ2U=0.95Table 2
The prescribed constraints and equivalent parameters for partly unacceptable planning (Vd and Dmean denote the percentage of volume receiving at least d [Gy] and mean dose, respectively).
Assigned region (colored region in Figure2(b))OrganConstraintEquivalent parameterPTV (red)PTVV75<5%(K1,b1U,ζ1U) with b1U=75 and ζ1U=0.95V70≥95%(K1,b1L,ζ1L) with b1L=70 and ζ1L=0.95OAR (blue)Spinal cordV30=0%(K2,b2U,ζ2U) with b2U=30 and ζ2U=1Dmean≤10GyK2,c2U with c2U=10OAR (green)ParotidV5<50%(K3,b3U,ζ3U) with b3U=5 and ζ3U=0.5Dmean≤10Gy(K3,c3U) with c3U=10


## 6. Experimental Results and Discussion

We used the iterative solutionszn, n=0,1,2,⋯, to the algorithms described by Equations (21) and (26) with h=1 and z0=x0 for solving the IMRT inverse problems of acceptable and remaining-one unacceptable plannings, respectively. The initial values of x0 were commonly chosen as xj0=0.1 for any j, unless otherwise specified.


### 6.1. Acceptable Planning

We first examined solving the problem of acceptable IMRT planning using the harder condition defined by AAPM, as shown in Table1, with the C-shape phantom in Figure 2(a). The AAPM report [35] says that “the goal of this problem set is probably not achievable and tests a system that is being pushed very hard,” but it is acceptable in the meaning of Definition 4. Actually, we obtained the dose-volume histograms (DVHs) showing that the dose distributions of PTV and OAR (core) perfectly satisfy all prescribed DVCs (see Figure 3) generated by using the solution z500. Figure 4 shows the iterative evolution of the divergence Vz defined in Equation (16), which is a Lyapunov function of the continuous analog described in Equation (13). The monotonic decrease of the divergence with solutions toward a fixed point in the set A confirms the theoretical result of Theorem 6 and guarantees an appropriate discretization of the differential equation.Figure 3
DVH of (a) PTV withK1,b1U,ζ1U and K1,b1L,ζ1L and (b) OAR with K2,b2U,ζ2U obtained for acceptable planning. The red right-angle corners indicate DVCs in Table 1.
(a)
K1,b1U,ζ1U, K1,b1L,ζ1L(b)
K2,b2U,ζ2UFigure 4
Time evolution of divergenceV for acceptable planning.The evolutions of the dose-volume proportion rates(30)I1UznI1,I1LznI1,I2UznI2,n=0,1,2,⋯,500,for DVCs with K1,b1U,ζ1U, K1,b1L,ζ1L, and K2,b2U,ζ2U are, respectively, shown in the upper panel of Figure 5. We see that each proportion rate behaves to become greater than the prescribed rate (red line in the graph) as the number of iteration steps increases. The dose distribution satisfying Equation (5) implies a partly dose-volume acceptability according to Definition 2 and is confirmed by each graph of binary functions R1Uzn, R1Lzn, and R2Uzn as shown in the lower panel of Figure 5. When all dose distributions become partly dose-volume acceptable, then every binary function is exactly zero after the iterative solution has reached a fixed point in A.Figure 5
Time evolutions of (a)∣I1U∣/I1 and R1U, (b) ∣I1L∣/I1 and R1L, and (c) ∣I2U∣/I2 and R2U for acceptable planning. A red line indicates the prescribed proportion rate.
(a)
K1,b1U,ζ1U(b)
K1,b1L,ζ1L(c)
K2,b2U,ζ2U


### 6.2. Remaining-One Unacceptable Planning

The proposed Procedure7 in Section 4 was experimentally demonstrated by applying it to a remaining-one unacceptable dose-volume and mean-dose constraint planning for the head and neck phantoms in Figure 2(b) with constraints given in Table 2.Figure6 shows DVHs obtained from a solution z1000 to the algorithm in Equation (26) with λ=1, N1=0, and z0=x0 for the totally unacceptable planning S defined in Equation (23). The resultant values V75=28% and V70=67% at PTV, V30=1% at OAR (spinal cord), and so on are not satisfied with the prescribed constraints except for Dmean=5.7 Gy at OAR (parotid). Because the achievement of constraints for PTV and OAR (spinal cord) are mandatory based on the overall objective, we tried to examine an inverse problem for optimizing the subsystem Se in Equation (25) using a fixed point z100∈Ae applied to the algorithm described in Equation (21). As shown in Figure 7 (which indicates DVHs), the dose distributions of PTV and OAR (spinal cord) fulfill the prescribed conditions as expected; however, the results of V5=100% and Dmean=50 Gy at OAR (parotid) are clinically unsatisfactory. For the perfect satisfaction of dose-volume and mean-dose constraints at PTV and OAR (spinal cord) being kept while the dose distribution at OAR (parotid) is reduced without causing the user’s trial and error, we performed an experiment according to the proposed procedure. The value of the parameter λ in Equation (21) used in Step 2 of Procedure 7 was set to 0.9, unless otherwise noted. The solution z1000 obtained after the procedure was used for calculating DVHs, as shown in Figure 8. We see that the dose distributed to OAR (parotid) was less than that shown in Figure 7(c). Actually, Dmean at OAR (parotid) can be decreased to 29 Gy from 50 Gy owing to the procedure. It is reported that a lower mean dose of the parotid gland is more effective for avoiding mouth dryness caused by salivary gland disorder [42]. Moreover, one can prevent a salivary disorder by reducing the mean doses to be applied to the parotid gland so that they are as small as possible [37, 38].Figure 6
DVH of (a) PTV withK1,b1U,ζ1U and K1,b1L,ζ1L, (b) OAR (spinal cord) with K2,b2U,ζ2U, and (c) OAR (parotid) with K3,b3U,ζ3U obtained for totally unacceptable planning. The red right-angle corners indicate DVCs in Table 2.
(a)
K1,b1U,ζ1U, K1,b1L,ζ1L(b)
K2,b2U,ζ2U(c)
K3,b3U,ζ3UFigure 7
DVH of (a) PTV withK1,b1U,ζ1U and K1,b1L,ζ1L, (b) OAR (spinal cord) with K2,b2U,ζ2U, and (c) OAR (parotid) with K3,b3U,ζ3U obtained for partly acceptable planning. The red right-angle corners indicate DVCs in Table 2.
(a)
K1,b1U,ζ1U, K1,b1L,ζ1L(b)
K2,b2U,ζ2U(c)
K3,b3U,ζ3UFigure 8
DVH of (a) PTV withK1,b1U,ζ1U and K1,b1L,ζ1L, (b) OAR (spinal cord) with K2,b2U,ζ2U, and (c) OAR (parotid) with K3,b3U,ζ3U obtained by the proposed procedure for remaining-one unacceptable planning. The red right-angle corners indicate DVCs in Table 2.
(a)
K1,b1U,ζ1U, K1,b1L,ζ1L(b)
K2,b2U,ζ2U(c)
K3,b3U,ζ3UFor evaluating the effectiveness of the proposed procedure quantitatively, we drew the graph in Figure9 showing the evolution of the distances
(31)Uez=R1Uz∑i1=1I1P1UK1zi1−K1zi12+R1Lz∑i1=1I1P1LK1zi1−K1zi12+R2Uz∑i2=1I2P2UK2zi2−K2zi22+S2Uz∑i2=1I2Q2UK2zi2−K2zi22,Urz=R3Uz∑i3=1I3P3UK3zi3−K3zi32+S3Uz∑i3=1I3Q3UK3zi3−K3zi32,each of which is defined as a squared L2 norm. The values of Uez and Urz become zero if and only if the solutions z to the subsystems Se and Sr are in the acceptable sets Ae and Ar, respectively. In the figure, the numbers of iteration n in the ranges 0,N1 and N1,1000 where N1=100 are, respectively, in Steps 1 and 2 of Procedure 7. We see that the iterative sequence znn=0N1 converges to the subspace Ae with decreasing Ueznn=0N1 by applying the process of Step 1 for obtaining an initial state of Step 2. In Step 2, Urznn=N1N1+27 is drastically decreased although the solution overflows from Ae, and Ueznn=N1N1+5 is slightly increased. The remaining iterations are for making the solution converge to Ae again. Note that the maximum of the sequence Urznn=N1+271000 results in a small value compared to the decrease caused in the preceding iterations.Figure 9
Time evolution of distancesUe (blue) and Ur (red) obtained using Procedure 7 for remaining-one unacceptable planning.We show how selecting the parameter value affects the performance with respect to decreasing the cost function. Figure10 illustrates the relation between λ and UrzN1+N2 where N2 denotes the iteration number at which the solution to Equation (26) with initial state zN1 has fallen into the partly acceptable set Ae. From the figure, we observe that a larger λ leads to a smaller value of the resultant Ur, while being accompanied with a larger number of iterations, N2, for the convergence.Figure 10
Relation betweenλ versus UrzN1+N2 (red dots) and N2 (green dots) on left and right axes, respectively.


## 6.1. Acceptable Planning

We first examined solving the problem of acceptable IMRT planning using the harder condition defined by AAPM, as shown in Table1, with the C-shape phantom in Figure 2(a). The AAPM report [35] says that “the goal of this problem set is probably not achievable and tests a system that is being pushed very hard,” but it is acceptable in the meaning of Definition 4. Actually, we obtained the dose-volume histograms (DVHs) showing that the dose distributions of PTV and OAR (core) perfectly satisfy all prescribed DVCs (see Figure 3) generated by using the solution z500. Figure 4 shows the iterative evolution of the divergence Vz defined in Equation (16), which is a Lyapunov function of the continuous analog described in Equation (13). The monotonic decrease of the divergence with solutions toward a fixed point in the set A confirms the theoretical result of Theorem 6 and guarantees an appropriate discretization of the differential equation.Figure 3
DVH of (a) PTV withK1,b1U,ζ1U and K1,b1L,ζ1L and (b) OAR with K2,b2U,ζ2U obtained for acceptable planning. The red right-angle corners indicate DVCs in Table 1.
(a)
K1,b1U,ζ1U, K1,b1L,ζ1L(b)
K2,b2U,ζ2UFigure 4
Time evolution of divergenceV for acceptable planning.The evolutions of the dose-volume proportion rates(30)I1UznI1,I1LznI1,I2UznI2,n=0,1,2,⋯,500,for DVCs with K1,b1U,ζ1U, K1,b1L,ζ1L, and K2,b2U,ζ2U are, respectively, shown in the upper panel of Figure 5. We see that each proportion rate behaves to become greater than the prescribed rate (red line in the graph) as the number of iteration steps increases. The dose distribution satisfying Equation (5) implies a partly dose-volume acceptability according to Definition 2 and is confirmed by each graph of binary functions R1Uzn, R1Lzn, and R2Uzn as shown in the lower panel of Figure 5. When all dose distributions become partly dose-volume acceptable, then every binary function is exactly zero after the iterative solution has reached a fixed point in A.Figure 5
Time evolutions of (a)∣I1U∣/I1 and R1U, (b) ∣I1L∣/I1 and R1L, and (c) ∣I2U∣/I2 and R2U for acceptable planning. A red line indicates the prescribed proportion rate.
(a)
K1,b1U,ζ1U(b)
K1,b1L,ζ1L(c)
K2,b2U,ζ2U


## 6.2. Remaining-One Unacceptable Planning

The proposed Procedure7 in Section 4 was experimentally demonstrated by applying it to a remaining-one unacceptable dose-volume and mean-dose constraint planning for the head and neck phantoms in Figure 2(b) with constraints given in Table 2.Figure6 shows DVHs obtained from a solution z1000 to the algorithm in Equation (26) with λ=1, N1=0, and z0=x0 for the totally unacceptable planning S defined in Equation (23). The resultant values V75=28% and V70=67% at PTV, V30=1% at OAR (spinal cord), and so on are not satisfied with the prescribed constraints except for Dmean=5.7 Gy at OAR (parotid). Because the achievement of constraints for PTV and OAR (spinal cord) are mandatory based on the overall objective, we tried to examine an inverse problem for optimizing the subsystem Se in Equation (25) using a fixed point z100∈Ae applied to the algorithm described in Equation (21). As shown in Figure 7 (which indicates DVHs), the dose distributions of PTV and OAR (spinal cord) fulfill the prescribed conditions as expected; however, the results of V5=100% and Dmean=50 Gy at OAR (parotid) are clinically unsatisfactory. For the perfect satisfaction of dose-volume and mean-dose constraints at PTV and OAR (spinal cord) being kept while the dose distribution at OAR (parotid) is reduced without causing the user’s trial and error, we performed an experiment according to the proposed procedure. The value of the parameter λ in Equation (21) used in Step 2 of Procedure 7 was set to 0.9, unless otherwise noted. The solution z1000 obtained after the procedure was used for calculating DVHs, as shown in Figure 8. We see that the dose distributed to OAR (parotid) was less than that shown in Figure 7(c). Actually, Dmean at OAR (parotid) can be decreased to 29 Gy from 50 Gy owing to the procedure. It is reported that a lower mean dose of the parotid gland is more effective for avoiding mouth dryness caused by salivary gland disorder [42]. Moreover, one can prevent a salivary disorder by reducing the mean doses to be applied to the parotid gland so that they are as small as possible [37, 38].Figure 6
DVH of (a) PTV withK1,b1U,ζ1U and K1,b1L,ζ1L, (b) OAR (spinal cord) with K2,b2U,ζ2U, and (c) OAR (parotid) with K3,b3U,ζ3U obtained for totally unacceptable planning. The red right-angle corners indicate DVCs in Table 2.
(a)
K1,b1U,ζ1U, K1,b1L,ζ1L(b)
K2,b2U,ζ2U(c)
K3,b3U,ζ3UFigure 7
DVH of (a) PTV withK1,b1U,ζ1U and K1,b1L,ζ1L, (b) OAR (spinal cord) with K2,b2U,ζ2U, and (c) OAR (parotid) with K3,b3U,ζ3U obtained for partly acceptable planning. The red right-angle corners indicate DVCs in Table 2.
(a)
K1,b1U,ζ1U, K1,b1L,ζ1L(b)
K2,b2U,ζ2U(c)
K3,b3U,ζ3UFigure 8
DVH of (a) PTV withK1,b1U,ζ1U and K1,b1L,ζ1L, (b) OAR (spinal cord) with K2,b2U,ζ2U, and (c) OAR (parotid) with K3,b3U,ζ3U obtained by the proposed procedure for remaining-one unacceptable planning. The red right-angle corners indicate DVCs in Table 2.
(a)
K1,b1U,ζ1U, K1,b1L,ζ1L(b)
K2,b2U,ζ2U(c)
K3,b3U,ζ3UFor evaluating the effectiveness of the proposed procedure quantitatively, we drew the graph in Figure9 showing the evolution of the distances
(31)Uez=R1Uz∑i1=1I1P1UK1zi1−K1zi12+R1Lz∑i1=1I1P1LK1zi1−K1zi12+R2Uz∑i2=1I2P2UK2zi2−K2zi22+S2Uz∑i2=1I2Q2UK2zi2−K2zi22,Urz=R3Uz∑i3=1I3P3UK3zi3−K3zi32+S3Uz∑i3=1I3Q3UK3zi3−K3zi32,each of which is defined as a squared L2 norm. The values of Uez and Urz become zero if and only if the solutions z to the subsystems Se and Sr are in the acceptable sets Ae and Ar, respectively. In the figure, the numbers of iteration n in the ranges 0,N1 and N1,1000 where N1=100 are, respectively, in Steps 1 and 2 of Procedure 7. We see that the iterative sequence znn=0N1 converges to the subspace Ae with decreasing Ueznn=0N1 by applying the process of Step 1 for obtaining an initial state of Step 2. In Step 2, Urznn=N1N1+27 is drastically decreased although the solution overflows from Ae, and Ueznn=N1N1+5 is slightly increased. The remaining iterations are for making the solution converge to Ae again. Note that the maximum of the sequence Urznn=N1+271000 results in a small value compared to the decrease caused in the preceding iterations.Figure 9
Time evolution of distancesUe (blue) and Ur (red) obtained using Procedure 7 for remaining-one unacceptable planning.We show how selecting the parameter value affects the performance with respect to decreasing the cost function. Figure10 illustrates the relation between λ and UrzN1+N2 where N2 denotes the iteration number at which the solution to Equation (26) with initial state zN1 has fallen into the partly acceptable set Ae. From the figure, we observe that a larger λ leads to a smaller value of the resultant Ur, while being accompanied with a larger number of iterations, N2, for the convergence.Figure 10
Relation betweenλ versus UrzN1+N2 (red dots) and N2 (green dots) on left and right axes, respectively.


## 7. Conclusion

For handling dose-volume and mean-dose constraints directly in the optimization of IMRT treatment inverse planning, we have proposed a novel iterative algorithm derived from the discretization of a continuous-time dynamical system. The proposed system has an advantage in that the stability of an equilibrium corresponding to the desired optimal solution to the inverse problem is able to be proven using the Lyapunov theorem. Through numerical experiments with the C-shape phantom (AAPM TG-119) for an acceptable planning and an extended C-shape phantom simulating the head and neck for remaining-one unacceptable planning, we confirmed that we can obtain an optimal solution and a nearly optimal solution located as close as possible to the remaining-one unacceptable constraint.Our approach presented in this paper enables us to develop an iterative method of IMRT optimization by discretizing a continuous-time system in which the global stability of a desired equilibrium is guaranteed based on the dynamical system theory. The advantage of the approach is due to the Lyapunov theorem that establishes the stability of equilibrium if there exists a Lyapunov function even for a hybrid dynamical system. The direct construction of iterative algorithms with a theoretical guarantee of convergence for a given objective function including inequalities with percentages is generally difficult. Although a drawback of the dynamical system approach is that finding a combination of a vector field and Lyapunov function is generally a hard problem, we succeeded in obtaining a proof for the acceptable set and designing an iterative optimization method for an IMRT treatment planning including dose-volume constraints.



---

*Source: 3096067-2020-09-03.xml*