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.
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 constraints directly 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 using the proposed optimizing system for the objective function such that the corresponding inverse problem is partly acceptable by excluding a function from the total objective function with an acceptable set with . Then, the next step is to move the state to decrease and 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 and , with an initial value in . The optimizing system can get a desired solution in that is close to the optimal subspace for the remaining-one constraint function .
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: indicates the th element of the vector ; and and indicate the th row and th 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, let 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 with 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 with indicating the closure of the open hypercube for some which represents the intensity in the corresponding beamlet, and which consists of all elements representing the dose deposited in the th voxel due to the unit intensity in the th beamlet; it is called the dose influence (or deposition) matrix. If there are volumes of PTV and OAR with and voxels, respectively, then includes and as subvectors, and includes and as submatrices. A similar definition can be applied to the existence of multiple PTVs and OARs.
Let and represent the lower and upper bounds on the dose delivered to PTV and OAR, respectively ( and ). Additionally, we define an upper dose bound on the dose delivered to PTV as to avoid excessively high doses inside the PTV. We also handle a mean dose denoted by as a real-valued function of the dose , which is defined by for . 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 and , respectively.
Definition 1. The IMRT planning system is consistent if the set is not empty; otherwise, it is inconsistent.
Definition 2. For the given and each set of dose volumes and their conditions , , and , where , , and 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 , , and , namely, each of the inequalities is satisfied for some , where indicates the number of elements in the set.
Definition 3. For the given and each set of doses and their conditions (, ) and (, ), 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 mean-dose constraints.
Definition 4. The IMRT planning system 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 if the system is acceptable.
For describing the system, let us define the functions , , , , and as follows:
for and . We consider the case where
with . To simplify the description, we also define
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 where for with . Where for and .
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 of the dynamical system in Equation (13), then the solution stays in for all .
Proof. Since the system can be written as with a function , we see that for any , the solution satisfies on the subspace where or . 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 of the system in Equation (13) with initial value being in for all .
Next, we prove the stability of equilibrium in the set , 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 in Equation (9) as an equilibrium set of Equation (13) is stable.
Proof. Any point is an equilibrium of Equation (13), so we can say that 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 . We obtain its derivative along the solution to the system in Equation (13). When is neither partly dose-volume nor mean-dose acceptable, the derivative is given by To treat any , 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 (, ) and (, ). The derivative is zero at . Thus, is a Lyapunov function with respect to . Consequently, the equilibrium set 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 for radiation beam weight for and with , where denotes a step size. When , , and all binary functions , , , , and are assumed to be one for simplicity, the th element of can be described as for . We see that the iterative formula in Equation (22) is similar to that of the expectation-maximization (EM) method in computed tomography. Namely, replacing with the measured projection and assuming and 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 .
4. 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 while satisfying the partly acceptable situation as in the subsystem . In the procedure, we use an initial value problem for the following difference equation of variables with initial state starting at : for , where the functions and are, respectively, defined by for . Where for and with a positive parameter . Note that one takes and then has , , for , when . We give the procedure as follows:
Procedure 7. For an IMRT planning system consisting of Equation (23), which is totally unacceptable and has a partly acceptable subsystem in Equation (25), the following shows a procedure for obtaining a solution located as close as possible to the remaining-one unacceptable constraint of in Equation (24) while satisfying the partly acceptable situation as in the subsystem .
Step 1. The first step of the procedure is to solve iterative solutions to the difference equation in Equation (21) for the subsystem and obtain a stable fixed point in the set illustrated in Figure 1. Obtaining 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 and an initial value and move the state , , while expecting the decrease of the objective function , as shown in Figure 1. Note that because the partly acceptable situation of the subsystem means that in Equation (26) at for , the iterative steps can forcibly decrease the sequence for at least a small value of . By introducing the function with a time-varying parameter in Equation (26), the iterative process has a time-dependent regularization effect making it possible to optimize multiple criteria. As the function tends to independently of along with time that has passed, the optimizing system can get a desired solution in that is close to the optimal subspace with respect to the objective function 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.” Figure 2 shows -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].
(a)
(b)
Table 1 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.
6. Experimental Results and Discussion
We used the iterative solutions , , to the algorithms described by Equations (21) and (26) with and for solving the IMRT inverse problems of acceptable and remaining-one unacceptable plannings, respectively. The initial values of were commonly chosen as for any , 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 . Figure 4 shows the iterative evolution of the divergence 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 confirms the theoretical result of Theorem 6 and guarantees an appropriate discretization of the differential equation.
(a) ,
(b)
The evolutions of the dose-volume proportion rates for DVCs with , , and 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 , , and 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)
(b)
(c)
6.2. Remaining-One Unacceptable Planning
The proposed Procedure 7 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.
Figure 6 shows DVHs obtained from a solution to the algorithm in Equation (26) with , , and for the totally unacceptable planning defined in Equation (23). The resultant values and at PTV, at OAR (spinal cord), and so on are not satisfied with the prescribed constraints except for 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 in Equation (25) using a fixed point 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 and 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 , unless otherwise noted. The solution 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, at OAR (parotid) can be decreased to Gy from 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].
(a) ,
(b)
(c)
(a) ,
(b)
(c)
(a) ,
(b)
(c)
For evaluating the effectiveness of the proposed procedure quantitatively, we drew the graph in Figure 9 showing the evolution of the distances each of which is defined as a squared norm. The values of and become zero if and only if the solutions to the subsystems and are in the acceptable sets and , respectively. In the figure, the numbers of iteration in the ranges and where are, respectively, in Steps 1 and 2 of Procedure 7. We see that the iterative sequence converges to the subspace with decreasing by applying the process of Step 1 for obtaining an initial state of Step 2. In Step 2, is drastically decreased although the solution overflows from , and is slightly increased. The remaining iterations are for making the solution converge to again. Note that the maximum of the sequence 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 where denotes the iteration number at which the solution to Equation (26) with initial state has fallen into the partly acceptable set . From the figure, we observe that a larger leads to a smaller value of the resultant , while being accompanied with a larger number of iterations, , for the convergence.
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.
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.
Acknowledgments
This research was partially supported by JSPS KAKENHI Grant Number 18K04169.