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

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.


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][4][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][7][8][9][10][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][15][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 constraints directly in IMRT treatment planning, which is an iterative algorithm derived from the discretization of a continuoustime dynamical system [22][23][24][25][26][27][28][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]. More-over, 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 A e using the proposed optimizing system for the objective function U e such that the corresponding inverse problem is partly acceptable by excluding a function U r from the total objective function with an acceptable set A r with A e ∩ A r = ϕ. Then, the next step is to move the state x to decrease U r ðxÞ and U e ðxÞ 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 U e and U r , with an initial value a e in A e . The optimizing system can get a desired solution in A e that is close to the optimal subspace A r for the remaining-one constraint function U r . This paper is organized as follows. In Section 2, 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 θ; Acceptable except for one constraint Remaining-one unacceptable constraint a e e r U e (x), U r (x) x ∈ Ω and ðΘÞ i and ðΘÞ j indicate the ith row and jth column vectors of the matrix Θ, respectively.

IMRT Planning
The radiation oncologist performs IMRT planning to determine the PTV to be irradiated and the OARs to be spared. With IMRT, let J 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 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 threedimensional 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 k ij 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 I 1 and I 2 voxels, respectively, then D includes D 1 ∈ R I 1 + and D 2 ∈ R I 2 + as subvectors, and K includes K 1 ∈ R I 1 ×J + and K 2 ∈ R I 2 ×J + as submatrices. A similar definition can be applied to the existence of multiple PTVs and OARs.
Let b L 1 and b U 2 represent the lower and upper bounds on the dose delivered to PTV and OAR, respectively (b L 1 > 0 and b U 2 ≥ 0). Additionally, we define an upper dose bound on the dose delivered to PTV as b U 1 > b L 1 to avoid excessively high doses inside the PTV. We also handle a mean dose denoted by kD m k as a real-valued function of the dose D m ∈ R I m + , which is defined by 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 c L 1 and c U 2 , respectively.
is not empty; otherwise, it is inconsistent.
Definition 2. For the given x ∈ Ω and each set of dose volumes and their conditions ðK are the prescribed proportion rates, the corresponding dose distribution is partly dose-volume acceptable if each number of elements of the index sets is, respectively, greater than the prescribed proportion of I 1 , I 1 , and I 2 , namely, each of the inequalities is satisfied for some x, where j·j indicates the number of elements in the set.
Definition 3. For the given x ∈ Ω and each set of doses and their conditions (K 1 , c L 1 ) and (K 2 , c U 2 ), the corresponding dose distribution is partly mean-dose acceptable if each of the inequalities is 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 meandose constraints.
Þg is acceptable if the following set is not empty: 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 variable x ∈ A if the system is acceptable. For describing the system, let us define the functions P U 1 ðD 1 Þ, P L 1 ðD 1 Þ, P U 2 ðD 2 Þ, Q L 1 ðD 1 Þ, and Q U 2 ðD 2 Þ as follows: for i 1 = 1, 2, ⋯, I 1 and i 2 = 1, 2, ⋯, I 2 . We consider the case where with I = 3I 1 + 2I 2 . To simplify the description, we also define

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 where for for i 1 = 1, 2, ⋯, I 1 and i 2 = 1, 2, ⋯, I 2 .
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 value x 0 ∈ Ω of the dynamical system in Equation (13), then the solution φðt, x 0 Þ stays in Ω for all t > 0.
Proof. Since the system can be written as dx j /dt = x j ð1 − γ −1 x j Þh j ðxÞ with a function h j , we see that for any j, the solution satisfies dφ j /dt ≡ 0 on the subspace where x j = 0 or x j = γ. 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, x 0 Þ of the system in Equation (13) with initial value x 0 ∈ Ω being in Ω for all t > 0.
Next, we prove the stability of equilibrium in the set A, 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, then A in Equation (9) as an equilibrium set of Equation (13) is stable.
Proof. Any point a ∈ 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 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 To treat any x ∈ Ω, define the index sets: Then, the calculation of the derivative is reduced to the following: The last inequality is supported by for partly mean-dose unacceptable (K 1 , c L 1 ) and (K 2 , c U 2 ). The derivative is zero at x = a ∈ Ω. Thus, VðxÞ 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 for j = 1, 2, ⋯, J and n = 0, 1, 2, ⋯ with zð0Þ = x 0 ∈ Ω, where h > 0 denotes a step size. When γ ⟶ ∞, h = 1, and all binary functions R U 1 , R L 1 , R U 1 , S 1 , and S U 2 are assumed to be one for simplicity, the jth element of z can be described as 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 PðKzðnÞÞ 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.

IMRT Optimization for Remaining-One Unacceptable Dose-Volume and Mean-Dose Constraint Planning
Consider an IMRT planning system where the subsystem 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 is acceptable. We show a procedure to get a solution located as close as possible to the remaining-one unacceptable constraint of S r while satisfying the partly acceptable situation as in the subsystem S e . In the procedure, we use an initial value problem for the following difference equation of variables zðnÞ with initial state zðN 1 Þ ∈ Ω starting at N 1 ≥ 0: for n = N 1 , N 1 + 1, ⋯, where the functions g j and δ j are, respectively, defined by for j = 1, 2, ⋯, J. Where for i 3 = 1, 2, ⋯, I 3 and with a positive parameter λ. Note that one takes λ n ⟶ 0 and then has g j ðzðnÞÞ ⟶ δ j , j = 1, 2, ⋯, J, for n ⟶ ∞, when λ < 1. We give the procedure as follows: Procedure 7. For an IMRT planning system S consisting of Equation (23), which is totally unacceptable and has a partly acceptable subsystem S e in Equation (25), the following shows a procedure for obtaining a solution located as close as possible to the remaining-one unacceptable constraint of S r in Equation (24) while satisfying the partly acceptable situation as in the subsystem S e .
Step 1. The first step of the procedure is to solve iterative solutions to the difference equation in Equation (21) for the subsystem S e and obtain a stable fixed point a e in the set A e illustrated in Figure 1. Obtaining a e ∈ A e 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 zðN 1 Þ = a e ∈ A e and move the state zðN 1 + nÞ, n = 1, 2, ⋯, N 2 , while expecting the decrease of the objective function U r ðzðN 1 + nÞÞ, as shown in Figure 1. Note that because the partly acceptable situation of the subsystem S e means that f j ðzÞ = σ j in Equation (26) at z = a e ∈ A e for j = 1, 2, ⋯, J, the iterative steps can forcibly decrease the sequence fU r ðzðnÞÞg N 1 +ν n=N 1 for at least a small value of ν. By introducing the function g j ðzðnÞÞ with a time-varying parameter λ n−N 1 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 zðnÞ along with time n that has passed, the optimizing system can get a desired solution in A e that is close to the optimal subspace A r with respect to the objective function U r for the remaining-one constraint.

Experimental Method
To evaluate the proposed method, we examine two IMRT treatment planning problems: "acceptable planning" and "remaining-one unacceptable planning." Figure 2 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][37][38][39][40][41][42]. Table 1 shows the DVCs for the C-shape phantom case, and we adopted a 9-field beam arrangement with every 40degree 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][4][5]. So we adopted a VMAT technique to perform complex treatment planning.

Experimental Results and Discussion
We used the iterative solutions zðnÞ, n = 0, 1, 2, ⋯, to the algorithms described by Equations (21) and (26) with h = 1 and zð0Þ = x 0 for solving the IMRT inverse problems of acceptable and remaining-one unacceptable plannings, respectively. The initial values of x 0 were commonly chosen as x 0 j = 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 Table 1, 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 zð500Þ. Figure 4 shows the iterative evolution of the divergence VðzÞ 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.
The evolutions of the dose-volume proportion rates for DVCs with ðK 1 , b U 1 , ζ U 1 Þ, ðK 1 , b L 1 , ζ L 1 Þ, and ðK 2 , b U 2 , ζ U 2 Þ 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 R U 1 ðzðnÞÞ, R L 1 ðzðnÞÞ, and R U 2 ðzðnÞÞ 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.
6.2. Remaining-One Unacceptable Planning. The proposed Procedure 7 in Section 4 was experimentally demonstrated by applying it to a remaining-one unacceptable dosevolume and mean-dose constraint planning for the head and neck phantoms in Figure 2(b) with constraints given in Table 2. Figure 6 shows DVHs obtained from a solution zð1000Þ to the algorithm in Equation (26) with λ = 1, N 1 = 0, and zð0Þ = x 0 for the totally unacceptable planning S defined in Equation (23). The resultant values V 75 = 28% and V 70 = 67% at PTV, V 30 = 1% at OAR (spinal cord), and so on are not satisfied with the prescribed constraints except for D mean = 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 S e in Equation (25) using a fixed point zð100Þ ∈ A e 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 V 5 = 100% and D mean = 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 zð1000Þ 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, D mean 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]. Assigned region (colored region in Figure 2(a)) Organ Constraint Equivalent parameter with b U 2 = 10 and ζ U 2 = 0:95 Assigned region (colored region in Figure 2 Figure 3: DVH of (a) PTV with ðK 1 , b U 1 , ζ U 1 Þ and ðK 1 , b L 1 , ζ L 1 Þ and (b) OAR with ðK 2 , b U 2 , ζ U 2 Þ obtained for acceptable planning. The red rightangle corners indicate DVCs in Table 1.    Figure 6: DVH of (a) PTV with ðK 1 , b U 1 , ζ U 1 Þ and ðK 1 , b L 1 , ζ L 1 Þ, (b) OAR (spinal cord) with ðK 2 , b U 2 , ζ U 2 Þ, and (c) OAR (parotid) with ðK 3 , b U 3 , ζ U 3 Þ obtained for totally unacceptable planning. The red right-angle corners indicate DVCs in Table 2.  Table 2. each of which is defined as a squared L 2 norm. The values of U e ðzÞ and U r ðzÞ become zero if and only if the solutions z to the subsystems S e and S r are in the acceptable sets A e and A r , respectively. In the figure, the numbers of iteration n in the ranges ½0, N 1 and ½N 1 , 1000 where N 1 = 100 are, respectively, in Steps 1 and 2 of Procedure 7. We see that the iterative sequence fzðnÞg N 1 n=0 converges to the subspace A e with decreasing fU e ðzðnÞÞg N 1 n=0 by applying the process of Step 1 for obtaining an initial state of Step 2. In Step 2, fU r ðzðnÞÞg N 1 +27 n=N 1 is drastically decreased although the solution overflows from A e , and fU e ðzðnÞÞg N 1 +5 n=N 1 is slightly increased. The remaining iterations are for making the solution converge to A e again. Note that the maximum of the sequence fU r ðzðnÞÞg 1000 n=N 1 +27 results in a small value compared to the decrease caused in the preceding iterations.
We show how selecting the parameter value affects the performance with respect to decreasing the cost function. Figure 10 illustrates the relation between λ and U r ðzðN 1 + N 2 ÞÞ where N 2 denotes the iteration number at which the solution to Equation (26) with initial state zðN 1 Þ has fallen into the partly acceptable set A e . From the figure, we observe  Table 2.  U e (z(n)), U r (z(n)) 3 Figure 9: Time evolution of distances U e (blue) and U r (red) obtained using Procedure 7 for remaining-one unacceptable planning. that a larger λ leads to a smaller value of the resultant U r , while being accompanied with a larger number of iterations, N 2 , for the convergence.

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.

Data Availability
All data used to support the findings of this study are included within the paper.

Conflicts of Interest
The authors declare no conflicts of interest regarding the publication of this paper.