Dose-volume constrained optimization in intensity-modulated radiation therapy treatment planning

Abstract

We present a novel optimization method to handle dose-volume constraints (DVCs) directly in intensity-modulated radiation therapy (IMRT) treatment planning based on the idea of continuous dynamical methods. Most of the conventional methods are constructed for solving inconsistent inverse problems with, e.g., dose-volume based objective functions, and one expects to obtain a feasible solution that minimally violates the DVCs. We introduce the concept of ‘acceptable’, meaning that there exists a nonempty set of radiation beam weights satisfying the given DVCs, and we resolve the issue that the objective and evaluation are different in the conventional planning approach. We apply the initial-value problem of the proposed dynamical system to an acceptable and inconsistent inverse problem and prove that the convergence to an equilibrium in the acceptable set of solutions is theoretically guaranteed by using the Lyapunov theorem. Indeed, we confirmed that we can obtain acceptable beam weights through numerical experiments using phantom data simulating a clinical setup for an acceptable and inconsistent IMRT planning system.

1 Introduction

Intensity-modulated radiation therapy (IMRT) [1, 2] is an effective radiotherapy technique that attacks cancer without inflicting serious damage to critical normal tissues by delivering radiation beams from many angles. For designing IMRT inverse plans, an optimization strategy is typically used to minimize an objective or cost function of radiation beam weights [3]. We widely utilize the objective function to achieve as high as possible proportions of the total volumes receiving doses greater than and less than the prescribed doses for planning target volumes (PTVs) and organs at risk (OARs), respectively. However, in practice, it is often impossible to fill all the volumes with permissible doses. Namely, the inverse problem generally becomes inconsistent, so the optimization problem has no optimum solution. Dose-volume constraints (DVCs) expressed as a percentage of the prescription dose are a natural way to specify the objective and have been the standard way of evaluating the treatment in practice. Most of the conventional methods are constructed for solving inconsistent inverse problems with, e.g., dose-based, dose-volume based, and biology-based objective functions [4–8], and one expects to obtain a feasible solution that minimally violates the DVCs. Therefore, planning for IMRT requires a trial-and-error process and human intervention [9] due to the difference between the optimization objective and the evaluation of planning results.

In this paper, we introduce the concept of ‘acceptable’ meaning that there exists a nonempty set of radiation beam weights satisfying the given DVCs on PTVs and OARs in the IMRT planning system and resolve the issue that the objective and evaluation are different in the conventional planning approach. That is, we present a novel optimization method to handle DVCs directly in IMRT treatment planning based on the idea of continuous dynamical methods [10–20]. We apply the initial-value problem of a dynamical system, which is described by nonlinear differential equations for an acceptable and inconsistent inverse problem, and prove that the convergence to an equilibrium [21] in the acceptable set of solutions is theoretically guaranteed by using the Lyapunov theorem [22]. Then we can obtain a set of acceptable beam weights over time with decreasing Kullback-Leibler divergence [23, 24] measure. We also show the fact that the intensities of all radiation beams are not negative and are less than a specific upper limit, which is achieved in accordance with a box-constrained optimization procedure. The proposed dynamical system extends our previously presented nonlinear continuous-time system [20] for solving a split feasibility problem [25–27]. A simulated dynamics approach [28] also uses a differential equation, but it is difficult to ascertain whether its solution has converged to the global optimum.

We examined numerical experiments using phantom data simulating a clinical setup for an acceptable and inconsistent IMRT planning system. The system proposed in this paper is compared with another dynamical system constructed for optimizing a consistent inverse problem. We found that the proposed system can provide an acceptable solution, but the other cannot.

2 IMRT planning

In this section, we give some definitions and notations. Assume that $x \in\bar{\Omega} \subset R^{J}$ is an unknown variable for radiation beam weight satisfying

$$ D = K x, $$

(1)

where $\bar{\Omega}$ indicates the closure of the open hypercube $\Omega = (0, \gamma)^{J}$ with $\gamma> 0$ applying an upper limit to all beam weights and $D \in R_{+}^{I}$ and $K \in R_{+}^{I \times J}$ denote the irradiated dose vector and the dose influence matrix, respectively, with $R_{+}$ being the set of nonnegative real numbers. If there are volumes of PTV and OAR with $I_{1}$ and $I_{2}$ voxels, respectively, then D includes $D_{1}\in R_{+}^{I_{1}}$ and $D_{2}\in R_{+}^{I_{2}}$ as subvectors, and K includes $K_{1}\in R_{+}^{I_{1}\times J}$ and $K_{2}\in R_{+}^{I_{2}\times J}$ as submatrices. A similar definition can be applied to the existence of multiple PTVs and OARs.

Let $D_{1}^{\mathrm{L}}$ and $D_{2}^{\mathrm{U}}$ represent the lower and upper bounds of doses for PTV and OAR, respectively ($D_{1}^{\mathrm{L}}> 0$ and $D_{2}^{\mathrm{U}}\ge0$). Additionally, we define an upper dose bound for PTV as $D_{1}^{\mathrm{U}}> D_{1}^{\mathrm{L}}$ to avoid excessively high dose values inside the PTV.

Definition 1

The IMRT planning system is consistent if the set

$$ {\mathcal{E}} = \bigl\{ e \in\Omega: D_{1}^{\mathrm {U}}\ge ( K_{1}e )_{i_{1}} \ge D_{1}^{\mathrm{L}}, \forall i_{1} \in \{1, \ldots, I_{1}\}, ( K_{2}e )_{i_{2}} \le D_{2}^{\mathrm{U}}, \forall i_{2} \in \{1, \ldots, I_{2}\} \bigr\} $$

(2)

is not empty; otherwise it is inconsistent.

Definition 2

For each set of dose volumes and their conditions $(\mathrm{PTV}, D_{1}^{\mathrm{U}}, \zeta_{1}^{\mathrm{U}})$, $(\mathrm{PTV}, D_{1}^{\mathrm{L}}, \zeta_{1}^{\mathrm{L}})$, and $(\mathrm{OAR}, D_{2}^{\mathrm{U}},\zeta _{2}^{\mathrm{U}})$, where $\zeta_{1}^{\mathrm{U}}\le1$, $\zeta_{1}^{\mathrm{L}}\le1$, and $\zeta_{2}^{\mathrm{U}}\le1$ are the prescribed proportion rates, the corresponding dose distribution is partly acceptable if there exists $x \in\Omega$ such that each number of elements of the index sets

$$\begin{aligned}& {\mathcal{I}}_{1}^{\mathrm{U}}(x) = \bigl\{ i_{1} \in \{1, \ldots, I_{1}\} : ( K_{1}x )_{i_{1}} \le D_{1}^{\mathrm{U}} \bigr\} , \\& {\mathcal{I}}_{1}^{\mathrm{L}}(x) = \bigl\{ i_{1} \in \{1, \ldots, I_{1}\} : ( K_{1}x )_{i_{1}} \ge D_{1}^{\mathrm{L}} \bigr\} ,\quad \mbox{and} \\& {\mathcal{I}}_{2}^{\mathrm {U}}(x) = \bigl\{ i_{2} \in \{1, \ldots, I_{2}\} : ( K_{2} x )_{i_{2}} \le D_{2}^{\mathrm{U}} \bigr\} \end{aligned}$$

(3)

is, respectively, greater than the prescribed proportion of $I_{1}$, $I_{1}$, and $I_{2}$, namely, each of the inequalities

$$\begin{aligned}& \#{\mathcal{I}}_{1}^{\mathrm{U}}(x) \ge \zeta _{1}^{\mathrm{U}} I_{1}, \\& \#{\mathcal{I}}_{1}^{\mathrm{L}}(x) \ge \zeta _{1}^{\mathrm{L}}I_{1},\quad \mbox{and} \\& \#{\mathcal{I}}_{2}^{\mathrm {U}}(x) \ge \zeta_{2}^{\mathrm{U}}I_{2} \end{aligned}$$

(4)

is satisfied for some x, where # indicates the number of elements in the set.

We define the IMRT planning system to be acceptable if there exists a common beam set such that dose distributions in PTVs and OARs are partly acceptable for all DVCs.

Definition 3

The IMRT planning system is acceptable if the following set is not empty:

$$ {\mathcal{A}} = \bigl\{ a \in \Omega: \#{\mathcal{I}}_{1}^{\mathrm{U}}(a) \ge\zeta_{1}^{\mathrm{U}}I_{1}, \#{\mathcal{I}}_{1}^{\mathrm{L}}(a) \ge\zeta_{1}^{\mathrm{L}}I_{1}, \# {\mathcal{I}}_{2}^{\mathrm{U}}(a) \ge\zeta_{2}^{\mathrm{U}}I_{2}\bigr\} . $$

(5)

Obviously, if the IMRT planning system is consistent, it is acceptable. We are interested in the situation where the system is inconsistent and acceptable. In this paper, the problem of dose-volume constrained optimization in IMRT planning is defined to obtain the unknown variable $x \in{\mathcal{A}}$ if the system is acceptable.

We define the projection P as follows:

$$ P : R_{+}^{I} \to R_{+}^{I};\qquad D = \left ( \begin{array}{@{}c@{}} D_{1} \\ D_{1} \\ D_{2} \end{array} \right ) \mapsto P(D), $$

(6)

where $I = 2 I_{1}+ I_{2}$ and

$$\begin{aligned}& \bigl(P(D)\bigr)_{i_{1}} = \left \{ \begin{array}{l@{\quad}l} ( D_{1} )_{i_{1}}, & \mbox{if } ( D_{1} )_{i_{1}} \le D_{1}^{\mathrm{U}}, \\ D_{1}^{\mathrm{U}}, & \mbox{otherwise}; \end{array} \right . \\& \bigl(P(D)\bigr)_{I_{1}+ i_{1}} = \left \{ \begin{array}{l@{\quad}l} ( D_{1} )_{i_{1}}, & \mbox{if } ( D_{1} )_{i_{1}} \ge D_{1}^{\mathrm{L}}, \\ D_{1}^{\mathrm{L}}, & \mbox{otherwise}; \end{array} \right . \\& \bigl(P(D)\bigr)_{2 I_{1}+ i_{2}} = \left \{ \begin{array}{l@{\quad}l} ( D_{2} )_{i_{2}}, & \mbox{if } ( D_{2} )_{i_{2}} \le D_{2}^{\mathrm{U}}, \\ D_{2}^{\mathrm{U}}, & \mbox{otherwise}; \end{array} \right . \end{aligned}$$

for $i_{1}= 1, 2, \ldots, I_{1}$ and $i_{2}= 1, 2, \ldots, I_{2}$.

We introduce the generalized Kullback-Leibler divergence [23] of two nonnegative vectors α and β:

$$ \operatorname{KL}(\alpha, \beta) = \sum_{\ell}\beta_{\ell}\log \frac{\beta_{\ell}}{\alpha_{\ell}} + \alpha_{\ell}- \beta_{\ell}, $$

(7)

where $\alpha_{\ell}$ and $\beta_{\ell}$ denote the ℓth elements of α and β, respectively. The divergence $\operatorname{KL}(\alpha, \beta)$ for the vectors α and β of nonnegative real numbers is nonnegative with $\operatorname{KL}(\alpha, \beta) = 0$ if and only if $\alpha= \beta$. This divergence is also called Csiszár’s I-divergence measure [24].

3 Dynamical system for dose-volume constrained optimization

This section provides an exact definition of the dose-volume constrained optimization method and its theoretical results. Consider an initial-value problem for the nonlinear differential equation

$$ \begin{aligned} &\frac{dx}{dt} = - X \bigl(U - \gamma^{-1} X\bigr) K^{\top}Q\bigl(\operatorname{Log}(K x) - \operatorname{Log}\bigl(P(K x)\bigr)\bigr), \\ &x(0) = x^{0}\in\Omega, \end{aligned} $$

(8)

where $X := \operatorname{diag}(x)$ indicates the diagonal matrix in which the diagonal entries starting in the upper left corner are the elements of x, U denotes the identity matrix, $\operatorname{Log}(y)$ is the vector-valued function $\operatorname{Log}(y) := (\log(y_{1}), \log(y_{2}), \ldots, \log(y_{I}))^{\top}$ of each element in vector $y = (y_{1}, y_{2}, \ldots, y_{I})^{\top}$, the superscript ⊤ stands for the transpose of a matrix or vector, $K \in R^{I\times J}$ is defined by

$$ K = \left ( \begin{array}{@{}c@{}} K_{1} \\ K_{1} \\ K_{2} \end{array} \right ), $$

and the projection Q is defined such that

$$ Q : R_{+}^{I} \to R_{+}^{I}; \qquad D = \left ( \begin{array}{@{}c@{}} D_{1} \\ D_{1} \\ D_{2} \end{array} \right ) \mapsto Q(D), $$

(9)

where $I = 2 I_{1}+ I_{2}$ and

$$\begin{aligned}& \bigl(Q(D)\bigr)_{i_{1}} = \left \{ \begin{array}{l@{\quad}l} 0, & \mbox{if } (\mathrm{PTV}, D_{1}^{\mathrm{U}}, \zeta_{1}^{\mathrm{U}}) \mbox{ is partly acceptable}, \\ ( D_{1} )_{i_{1}}, & \mbox{otherwise}; \end{array} \right . \\& \bigl(Q(D)\bigr)_{I_{1}+ i_{1}} = \left \{ \begin{array}{l@{\quad}l} 0, & \mbox{if }(\mathrm{PTV}, D_{1}^{\mathrm{L}}, \zeta_{1}^{\mathrm{L}}) \mbox{ is partly acceptable}, \\ ( D_{1} )_{i_{1}}, & \mbox{otherwise}; \end{array} \right . \\& \bigl(Q(D)\bigr)_{2 I_{1}+ i_{2}} = \left \{ \begin{array}{l@{\quad}l} 0, & \mbox{if }(\mathrm{OAR}, D_{2}^{\mathrm{U}}, \zeta_{2}^{\mathrm{U}})\mbox{ is partly acceptable}, \\ ( D_{2} )_{i_{2}}, & \mbox{otherwise}; \end{array} \right . \end{aligned}$$

for $i_{1}= 1, 2, \ldots, I_{1}$ and $i_{2}= 1, 2, \ldots, I_{2}$.

We give theoretical results for the behavior of the solution to the dynamical system in Eq. (8). First, we show that all solutions stay inside the hypercube.

Proposition 1

If we choose the initial value $x^{0} \in\Omega$ in the dynamical system in Eq. (8), then the solution $\varphi(t, x^{0})$ stays in Ω for all $t > 0$.

Proof

Since the system can be written as $dx_{j}/dt = - x_{j} (1 - \gamma^{-1} x_{j}) (K^{\top})_{j} Q(\operatorname {Log}(K x) - \operatorname{Log}(P(K x)))$, we see that for any j the solution satisfies $d\varphi_{j}/dt \equiv0$ on the subspace where $x_{j} = 0$ or $x_{j} = \gamma$. 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 $\varphi(t, x^{0})$ of the system in Eq. (8) with initial value $x^{0} \in\Omega$ being in Ω for all $t > 0$. □

Next, we prove the stability of an equilibrium in the set ${\mathcal {A}}$, which corresponds to the desired radiation beam weights. Namely, the existence of a Lyapunov function for the system in Eq. (8) guarantees the stability of the equilibrium set [21].

Theorem 1

If the IMRT planning system is acceptable, then the equilibrium set of Eq. (8) is stable.

Proof

Any point $a \in{\mathcal{A}}$ in Eq. (5) is an equilibrium of Eq. (8), which means that ${\mathcal{A}}$ is an equilibrium set as a union of equilibria. Consider a Lyapunov candidate function defined in the set $\bar{\Omega}$ as

$$\begin{aligned} W(x) =& \operatorname{KL}(Ka, Kx) \\ =& \sum_{i=1}^{I} (Kx)_{i} \log \frac{(Kx)_{i}}{(Ka)_{i}} + (Ka)_{i} - (Kx)_{i} \\ =& \sum_{i=1}^{I} \int _{(Ka)_{i}}^{(Kx)_{i}} \log\frac{y}{(Ka)_{i}}\, dy, \end{aligned}$$

which is positive definite with respect to the point $a \in{\mathcal{A}}$.

Note that there exists a diagonal matrix $\tilde{X} \in R_{+}^{J\times J}$ such that

$$ {\tilde{X}}^{2} = X \bigl(U - \gamma^{-1}X\bigr), $$

(10)

where the diagonal elements are positive real numbers.

Then we obtain the derivative of W along the solution $\varphi(t, x^{0})$ with an arbitrary initial state $x^{0} \in \Omega$ of the system in Eq. (8), which can be written as $dx/dt = f(x)$, in the following expression:

$$\begin{aligned} {\biggl.\frac{d}{dt}{W}( \varphi)\biggr|_{\text{(8)}} } =&\frac{\partial W}{\partial x} f(\varphi) \\ =& -\bigl(\operatorname{Log}(K\varphi) - \operatorname{Log}(Ka) \bigr)^{\top}K \Phi\bigl(U - \gamma^{-1} \Phi\bigr) K^{\top} Q\bigl(\operatorname{Log}(K\varphi) - \operatorname{Log} \bigl(P(K \varphi)\bigr) \bigr) \\ =& -\bigl(\operatorname{Log}(K \varphi) - \operatorname{Log}(Ka) \bigr)^{\top}K \tilde{\Phi} (K \tilde{\Phi})^{\top}Q\bigl( \operatorname{Log}(K\varphi) - \operatorname{Log}\bigl(P(K\varphi)\bigr) \bigr) \\ =& -\bigl((K \tilde{\Phi})^{\top}\bigl(\operatorname{Log}(K \varphi) - \operatorname{Log}(Ka)\bigr) \bigr)^{\top}\bigl((K \tilde{ \Phi})^{\top}Q\bigl(\operatorname{Log}(K \varphi) -\operatorname{Log} \bigl(P(K\varphi)\bigr)\bigr)\bigr) \\ =& -\sum_{i_{1}\in \bar{\mathcal{I}}_{1}^{\mathrm{U}}( \varphi)} \bigl\lVert (K_{1}\tilde{ \Phi})_{i_{1}} \bigr\rVert _{2}^{2} \cdot \bigl( \log \bigl((K_{1} \varphi)_{i_{1}}\bigr) - \log\bigl((K_{1}a)_{i_{1}} \bigr) \bigr) \\ &{}\cdot \bigl( \log\bigl((K_{1}\varphi)_{i_{1}} \bigr) - \log\bigl(D_{1}^{\mathrm{U}}\bigr) \bigr) \\ &{} -\sum_{i_{1}\in\bar{\mathcal{I}}_{1}^{\mathrm{L}}(\varphi)} \bigl\lVert (K_{1}\tilde{ \Phi})_{i_{1}} \bigr\rVert _{2}^{2} \cdot \bigl( \log \bigl((K_{1}\varphi)_{i_{1}}\bigr) - \log\bigl((K_{1}a)_{i_{1}} \bigr) \bigr) \\ &{}\cdot \bigl( \log\bigl((K_{1} \varphi)_{i_{1}} \bigr) -\log\bigl(D_{1}^{\mathrm{L}}\bigr) \bigr) \\ &{} -\sum_{i_{2}\in \bar{\mathcal{I}}_{2}^{\mathrm{U}}( \varphi)} \bigl\lVert (K_{2} \tilde{ \Phi})_{i_{2}} \bigr\rVert _{2}^{2} \cdot \bigl( \log \bigl((K_{2} \varphi)_{i_{2}}\bigr) - \log\bigl((K_{2}a)_{i_{2}} \bigr) \bigr) \\ &{}\cdot \bigl( \log\bigl((K_{2}\varphi)_{i_{2}}\bigr) - \log\bigl(D_{2}^{\mathrm{U}}\bigr) \bigr) \\ \le& -\sum_{i_{1}\in\bar{\mathcal{I}}_{1}^{\mathrm{U}}(\varphi)} \bigl\lVert (K_{1} \tilde{ \Phi})_{i_{1}} \bigr\rVert _{2}^{2} \cdot \bigl( \log\bigl((K_{1} \varphi)_{i_{1}}\bigr) - \log \bigl(D_{1}^{\mathrm{U}}\bigr) \bigr)^{2} \\ &{} -\sum_{i_{1}\in \bar{\mathcal{I}}_{1}^{\mathrm{L}}( \varphi)} \bigl\lVert (K_{1} \tilde{\Phi})_{i_{1}} \bigr\rVert _{2}^{2} \cdot \bigl( \log\bigl((K_{1} \varphi)_{i_{1}}\bigr) - \log \bigl(D_{1}^{\mathrm{L}}\bigr) \bigr)^{2} \\ &{} -\sum_{i_{2}\in\bar{\mathcal{I}}_{2}^{\mathrm{U}}( \varphi)} \bigl\lVert (K_{2} \tilde{ \Phi})_{i_{2}} \bigr\rVert _{2}^{2} \cdot \bigl(\log \bigl((K_{2} \varphi)_{i_{2}}\bigr) - \log\bigl(D_{2}^{\mathrm{U}} \bigr) \bigr)^{2} \\ =& - \bigl\lVert (K \tilde{ \Phi})^{\top}\bigl(\operatorname{Log}(K \varphi) - \operatorname{Log}\bigl(P(K\varphi)\bigr) \bigr) \bigr\rVert _{2}^{2} \\ \le& 0. \end{aligned}$$

Here Φ denotes $\operatorname{diag}(\varphi)$, $\tilde{\Phi}$ satisfies ${\tilde{\Phi}}^{2} = \Phi(U - \gamma^{-1}\Phi)$, and

$$\begin{aligned}& \bar{\mathcal{I}}_{1}^{\mathrm{U}}(x) = \left \{ \begin{array}{l@{\quad}l} \phi, & \mbox{if } (\mathrm{PTV}, D_{1}^{\mathrm{U}}, \zeta_{1}^{\mathrm{U}}) \mbox{ is partly acceptable}, \\ \{1, \ldots, I_{1}\} \setminus{\mathcal{I}}_{1}^{\mathrm{U}}(x), & \mbox{otherwise}; \end{array} \right . \\& \bar{\mathcal{I}}_{1}^{\mathrm{L}}(x) = \left \{ \begin{array}{l@{\quad}l} \phi, & \mbox{if }(\mathrm{PTV}, D_{1}^{\mathrm{L}}, \zeta_{1}^{\mathrm{L}})\mbox{ is partly acceptable}, \\ \{1, \ldots, I_{1}\} \setminus{\mathcal{I}}_{1}^{\mathrm{L}}(x), & \mbox{otherwise}; \end{array} \right . \\& \bar{\mathcal{I}}_{2}^{\mathrm{U}}(x) =\left \{ \begin{array}{l@{\quad}l} \phi, & \mbox{if }(\mathrm{OAR},D_{2}^{\mathrm{U}}, \zeta_{2}^{\mathrm{U}})\mbox{ is partly acceptable}, \\ \{1, \ldots, I_{2}\} \setminus{\mathcal{I}}_{2}^{\mathrm{U}}(x), & \mbox{otherwise}; \end{array} \right . \end{aligned}$$

with

$$ \bar{\mathcal{I}}_{1}^{\mathrm{U}}(x) \cap\bar{\mathcal {I}}_{1}^{\mathrm{L}}(x) = \phi, \quad \forall x \in\Omega. $$

The derivative is zero at $x = a \in \bar{\Omega}$. Thus, $W(x)$ is a Lyapunov function with respect to a. Consequently, the equilibrium set ${\mathcal{A}}$ is stable. □

The theoretical results show that an element in the acceptable set can be obtained by applying the initial-value problem of the hybrid dynamical system with piecewise continuous vector fields in Eq. (8). To be acceptable, the vector of beams or the solution $\varphi(t, x^{0})$ must behave appropriately for the volume percentage in the dose-volume constraints, roughly speaking, in the following manner. Let us assume, e.g., the dose $D_{1}= K_{1} \varphi(t, x^{0})$ for PTV does not satisfy the upper bound constraint $(D_{1}^{\mathrm{U}},\zeta_{1}^{\mathrm{U}})$, i.e., it is not partly acceptable. The solution $\varphi(t, x^{0})$ causes the dose to change its distributions along the gradient of the Kullback-Leibler divergence measure so that $D_{1}$ can satisfy the constraint. The same manner is applied to the lower and upper bounds of doses for PTV and OAR, respectively. When all dose distributions in PTVs and OARs become partly acceptable, the vector field is entirely zero after the solution has converged to an equilibrium.

4 Materials for numerical experiments

We examine the problem of IMRT treatment planning with DVCs for a $512 \times512$ image as shown in Figure 1. This phantom image is a slice of computed tomography (CT) images obtained from the Visible Human Project [29]. We assume that the prostate contains a cancerous malignant tumor and the other organs are normal. To treat the prostate cancer, we designated the blue region containing the prostate as a PTV and assigned the regions colored green, red, cyan, and magenta to OARs (OAR1 to OAR4). The numbers of voxels in PTV, OAR1, OAR2, OAR3, and OAR4 are 1,093, 1,971, 2,061, 818, and 2,476, respectively, and then we have $I = 9\text{,}512$.

Figure 1

CT image and designated PTV and OARs regions. (a) $512 \times512$ CT image of male lower abdomen and (b) colored regions designated as PTV and OARs. Blue region is PTV, and green, red, cyan, and magenta regions are OARs.

Table 1 shows the DVCs that were determined by referencing the ultra high prescribed dose (86.4 Gy) IMRT case for PTV [30, 31] and the normal prescribed dose (74-78 Gy) IMRT case for OARs [32, 33]. While the ultra high dose IMRT is very efficient to treat the prostate cancer [31], the DVCs to OARs used in normal IMRT cases cannot be applied ordinarily due to the increase of irradiated OARs with escalating dose. Despite this, in our experiments, all upper dose limits for the OARs were decreased by subtracting 20-35 Gy from the typical dose values in normal IMRT cases for creating a higher quality IMRT plan. In fact, the DVCs shown in Table 1 are very stringent.

Table 1 Dose-volume constraint or equivalent parameters for partly acceptable dose distribution in our example

For the phantom imitating prostate cancer, in accordance with the Memorial Sloan Kettering Cancer Center (MSKCC) protocol of the report [30, 31], we applied five-field irradiation at angles of 45, 135, 195, 270, and 345 degrees; these angles were measured counterclockwise from the horizontal line that passes through the center of the CT image in Figure 1. Since we should find that the intensity of radiation beams affects only the PTV and OARs, we obtain $J = 806$. In the calculation of K, we did not take scattering of radiation [34] into consideration to simplify the dose calculation.

The IMRT planning derived from the above mentioned construction is inconsistent and acceptable in the meaning of the definitions in Section 2.

5 Experimental results

For the IMRT planning system that is not consistent but is acceptable, we compared the proposed system of Eq. (8) with a continuous-time dynamical system for optimizing a consistent inverse problem, which is described by

$$ \begin{aligned} &\frac{dx}{dt} = - X \bigl(U - \gamma^{-1} X\bigr) K^{\top}\bigl(\operatorname{Log}(K x) - \operatorname{Log}\bigl(P(K x)\bigr)\bigr), \\ &x(0) = x^{0} \in\Omega. \end{aligned} $$

(11)

Note that, if the IMRT planning system is consistent, ℰ in Eq. (2) as an equilibrium set of Eq. (11) is stable. This fact can be proved in a similar manner as Theorem 1.

As with conventional IMRT planning optimization, the method of the split feasibility problem as well as the gradient methods such as Newton’s method and the conjugate gradient method are known to be used for objective function minimization. All of these methods for searching for feasibility solutions are based on iterative procedure. However, the differential equation described in Eq. (11) is a continuous analog of the CQ-algorithm [26], which is designed for implementing the split feasibility problem with the advantage of not calculating a matrix inversion or Hessian over the conventional gradient methods, and is suitable for comparing and verifying the mechanism of obtaining acceptable solutions to the proposed continuous system in Eq. (8).

The common parameters and conditions for numerical simulation are as follows. For numerical integration of differential equations, we used the solver ode15s provided by MATLAB (MathWorks, Natick, USA). The parameter γ was set to 100, and the initial values of the state variables at $t = 0$ were chosen as $x_{j}^{0} = 10$ for any j. The upper beam limit γ is defined to avoid machine restrictions, and we are able to confirm that the IMRT planning system cannot reach any acceptable solution if the second inequality of Eq. (4) with x replaced by the vector containing each element of γ is not satisfied.

We first tried to solve the IMRT planning problem by using the continuous-time dynamical system (11). Figure 2 shows the dose-volume histograms (DVHs) generated from the solution at $t = 2\text{,}900$ to Eq. (11). Each red right-angle corner in the figures indicates the DVC shown in Table 1. These DVHs show that the dose distributions do not satisfy three DVCs, $(\mathrm{PTV}, 82.1, 0.87)$, $(\mathrm{OAR}1, 30, 0.9)$, and $(\mathrm{OAR}4, 30, 0.9)$ because the DVH on PTV passes the lower-left side of the upper red right-angle corner and the DVHs on OAR1 and OAR4 pass the right side of each red right-angle corner, i.e., a failed IMRT plan was obtained. Furthermore, when we tried several different initial values $x^{0}$, we did not obtain an IMRT plan such that all the given DVCs are satisfied. The results demonstrated that this IMRT planning problem is quite a challenge to solve.

Figure 2

DVH obtained by method using dynamical system in Eq. ( 11 ). (a) PTV, (b) OAR1, (c) OAR2, (d) OAR3, and (e) OAR4. Red right-angle corner indicates DVC.

On the other hand, in spite of the difficult conditions, the proposed system (8) can easily produce an IMRT plan where all dose distributions are acceptable. Figure 3 shows the DVHs generated from the solution at $t = 3.5$ to Eq. (8). We can confirm that all volumes of PTV and OARs fulfill the given DVCs. Namely, the solution converges to an equilibrium in set ${\mathcal{A}}$. In particular, compared with the DVHs in Figures 2(c) and 2(d), the spreads between the DVHs and the red right-angle corners in Figures 3(c) and 3(d) were decreased; therefore, adequate doses to the PTV were realized as shown in Figure 3(a).

Figure 3

DVH obtained by proposed method using dynamical system in Eq. ( 8 ). (a) PTV, (b) OAR1, (c) OAR2, (d) OAR3, and (e) OAR4. Red right-angle corner indicates DVC.

To illustrate the mechanism of the behavior of solutions toward an equilibrium in the acceptable set, the time courses showing the variation of Q in Eq. (9) for some DVCs are drawn in Figure 4. The values of 0 and 1 indicate the instantaneous states where the dose distributions are partly acceptable and not acceptable, respectively, for the corresponding DVCs. As seen from Figures 4(d) and 4(e), it is interesting that graphs for $(\mathrm{OAR}2, 20, 0.4)$ and $(\mathrm{OAR}3, 20, 0.6)$ transiently oscillate between 0 and 1 and tend to zero as a steady state. When $(\mathrm{PTV}, 95.9, 0.98)$ is partly acceptable at $t = 3.385$ as shown in Figure 4(b), all dose distributions become acceptable.

Figure 4

Time course showing variation of Q . (a) $(\mathrm{PTV}, 82.1, 0.87)$, (b) $(\mathrm{PTV}, 95.9, 0.98)$, (c) $(\mathrm{OAR}1, 30, 0.9)$, (d) $(\mathrm{OAR}2, 20, 0.4)$, (e) $(\mathrm{OAR}3, 20, 0.6)$, and (f) $(\mathrm{OAR4}, 30, 0.9)$. Values of 0 and 1 indicate states of being partly acceptable and not acceptable, respectively.

6 Discussion and conclusion

The purpose of this paper is to present a dose-volume constrained optimization method of IMRT treatment planning using the initial-value problem of the continuous-time dynamical system. First, let us discuss the theoretical results obtained by using the proposed method. Theorem 1 guarantees that the value of the Lyapunov function W, which corresponds to an objective or cost function, monotonically decreases along a solution obtained by the proposed method when an IMRT planning system is acceptable, which is mathematically defined for the existence of radiation beams satisfying the given DVCs. Additionally, when we choose an initial state in Ω, our method is able to find a box-constrained solution within Ω that is supported by Proposition 1.

Then we discuss the convergence to an equilibrium in the acceptable set of solutions in numerical experiments. The behavior can be explained qualitatively as follows. When the dose distribution for a PTV or an OAR becomes partly acceptable as time passes, the value of the corresponding function Q and the term added to the vector field (the right-hand side of Eq. (8)), which is required to make it acceptable, are zero in the dynamical system. Therefore, it is expected that the dynamics will make it easier for other OARs and PTVs to work to be partly acceptable. The function Q plays an essential role to control the cooperation among terms for DVCs in the vector field. If all dose distributions in PTVs and OARs become partly acceptable, then the vector field is entirely zero after the solution has converged to an equilibrium corresponding to the desired radiation beam weights.

In accordance with the theoretical result, we can take any initial state $x^{0}$ in Ω. Indeed, when performing numerical experiments with several other initial values $x_{j}^{0}$ for all j in the range $(0, \gamma)$, similar results of convergence to equilibria in the acceptable set ${\mathcal{A}}$ were obtained for all trials. Note that our goal is to reach an arbitrary element in ${\mathcal{A}}$, although an optimal solution is, in general, an isolated point in the state space in the usual definition of optimization problems. Nevertheless, it is significant to consider an optimization problem restricted to the set ${\mathcal{A}}$ from the viewpoint of another object. Regarding this, the following property of solutions was observed experimentally: the mean value of the state variables in an acceptable equilibrium after obtaining convergence was relatively small when the initial state $x_{j}^{0}$ was taken as a small value. In practical IMRT planning, the radiation beam weights are expected to be as weak as possible for reducing normal tissue doses. Hence, we can say that we should choose a small initial value so long as the selection of the initial state in Ω results in the same convergence to an element in ${\mathcal{A}}$.

A discretization of the differential equation in Eq. (11) using the Euler method leads to the iterated CQ-algorithm; therefore, the continuous system has the same convergence property on solutions as that of the difference equation. That is, the continuous system possesses a good characteristic that its solutions converge asymptotically to an equilibrium if the given inverse problem is consistent; otherwise, a close solution to the complete solution can be found [35]. By extending the continuous system, we have presented a hybrid dynamical system that can change the vector fields depending on partly acceptable conditions. For the inverse problem dealing with the discontinuous constraints defined by proportion rates, the global stability of solutions of the hybrid system can be proved theoretically using the tool of the Lyapunov theorem.

Although the proposed method of using dynamical systems with piecewise continuous vector fields is well designed for dose-volume constrained optimization, it requires numerical integration with a high computational cost to solve solutions to differential equations. In future work, we will attempt to construct an iterative method formulated by discretizing our differential equations for reducing the costs.