Dose calculations

During radiotherapy, the patient is irradiated with a predefined beam set denoted by B. Each beam in this study consists of m rows and n columns of beamlets, with each beamlet size being 1 x 1 cm. The set of all generated deliverable apertures is denoted by K, and the weight of the aperture κ is y_κ. The beamlets in set A_κ are delivered to the patient via aperture κ, and the treatment involves S structures, where each structure s ($s=\mathrm{1,}\dots \mathrm{,}S$) comprises vs voxels. The element in the deposition matrix is the deposition coefficient Wijs, which represents the dose received by j ($j=\mathrm{1,}\dots \mathrm{,}{v}_s$) in structure s from the beamlet i ($i\in A_{\kappa }$) allowed to pass through the aperture κ at unit intensity. The dose D_js can be expressed as follows: $D_{js}={\displaystyle \sum _{\kappa \in K}\left({\displaystyle \sum _{i\in A_{\kappa }}{W}_{ijs}}\right)}{y}_{\kappa }\mathrm{.}$ (1)

Column generation

In this study, the optimization problem is constructed as $\min F\left(D\right)=\min {\displaystyle \sum _{s=1}^S{\displaystyle \sum _{\xi =1}^{N_s}{F}_{\xi s}}}\left(D_s\right)\mathrm{,}$ (2) $\text{s}\text{.t}\mathrm{.}{\displaystyle \sum _{\kappa \in K}\left({\displaystyle \sum _{i\in A_{\kappa }}{W}_{ijs}}\right)}{y}_{\kappa }=D_{js}\mathrm{,}$ (3) $y_{\kappa }\ge \mathrm{0,}\kappa \in K\mathrm{.}$ (4) F(D) is the objective function of the optimization problem. For s, Ns sub-objective functions $F_{\xi s}\left(D_s\right)$ are used to constrain the received dose Ds. An excellent description of column generation was provided by Romeijn et al. [5]. In one iteration, a new aperture shape is generated by solving the pricing problem and is accepted into the apertures set. The weights of the apertures set are optimized as the master problem. The limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm for bound-constrained optimization [12-14] was employed to optimize the weights of the apertures generated in this study. The optimization was terminated when either the treatment plan met the requirements of the planner or the iteration number reached its limit.

When an aperture is generated, its cost should be calculated [5]. The pricing problem is $\underset{\kappa \in K}{{\displaystyle \min }}{\displaystyle \sum _{i\in A_{\kappa }}\left({\displaystyle \sum _{s=1}^S{\displaystyle \sum _{j=1}^{v_s}{W}_{ijs}{\pi }_{js}}}\right)}\mathrm{,}$ (5) where ${\pi }_{js}=\frac{\partial F_{js}\left(D_{js}\right)}{\partial D_{js}}$ is the Lagrange multiplier under the Karush–Kuhn–Tucker (KKT) condition. The interdigitation of the MLC was not allowed for the proposed optimization method, and the network flow could be employed to solve the pricing problem with this constraint. For beam $l\in B$, in the r th ($r=\mathrm{1,}\dots \mathrm{,}m$) row of the beamlets, $c_1\left(c_1=\mathrm{0,}\dots \mathrm{,}n\right)$ is used to mark the last beamlet of row r covered by the left leaf, and $c_2\left(c_2=\mathrm{1,}\dots \mathrm{,}n+1\right)$ is used to mark the first beamlet of row r covered by the right leaf. According to Eq. (5), the cost of $\left(r\mathrm{,}{c}_1\mathrm{,}{c}_2\right)$ is the sum of the gradient elements not covered by the leaves in the gradient map.

Therefore, the radiation beam was decomposed into rectangular beamlets in this study. Each beamlet gradient was calculated and the gradient of the corresponding beamlet was arranged according to the position of the beamlet in the beam, forming the aperture gradient map. An aperture shape that can improve the objective function was generated according to the gradient map. This generated aperture shape was equivalent to the negative-gradient descent direction for the optimization process. The search was performed along the negative-gradient descent direction, which rapidly reduced the function value. However, this did not indicate that the convergence speed of the steepest descent method was high. The sawtooth effect implies that the search direction of a negative gradient is not necessarily the fastest descent direction in the global range.

Conjugate gradient modulation

Compared with Newton’s and quasi-Newton methods, the conjugate gradient method uses simpler calculations, and its convergence speed is higher than that of the steepest descent method. To speed up the optimization process and improve the optimization quality of column generation, the conjugate gradient method was used to modulate the gradient to generate the aperture shape.

2.3.2

Search direction based on conjugate gradient modulation

The decrease in the objective function value of the column generation methods based on conjugate gradient directions with different coefficients in Eq. (7) was observed using the same objective function in a case of head and neck cancer (Fig. 1).

Fig. 1Full-screen View

Decrease in the objective function value of the above six column generation methods based on the conjugate gradient direction.

The objective function value of the column generation method based on the PRP conjugate gradient direction initially exhibited the fastest decline. During the second half of the iterative process, the objective function value of the column generation method based on the HS conjugate gradient direction exhibited the fastest convergence. Based on the above observations, an aperture shape optimization method, PRP-HS, based on modulation of the PRP and HS conjugate gradients, was proposed.

Two conjugate gradient descent directions were used to determine the search direction of the aperture shape. The PRP conjugate gradient direction had priority for decision-making, which decreased as the optimization proceeded, and the HS conjugate gradient direction increased in decision-making priority. The expression for $d_k^{\text{PRP}-\text{HS}}$ is $d_k^{\text{PRP}-\text{HS}}=\frac{1}{k}{d}_k^{\text{PRP}}+\left(1-\frac{1}{k}\right)d_k^{\text{HS}}\mathrm{.}$ (8) This formula can be rewritten as $\begin{array}{l}{d}_k^{\text{PRP}-\text{HS}}\\ =\frac{1}{k}\left(-g_k+{\beta }_k^{\text{PRP}}{d}_{k-1}^{\text{PRP}}\right)+\left(1-\frac{1}{k}\right)\left(-g_k+{\beta }_k^{\text{HS}}{d}_{k-1}^{\text{HS}}\right)\\ =-g_k+\left(\frac{1}{k}{\beta }_k^{\text{PRP}}{d}_{k-1}^{\text{PRP}}+\left(1-\frac{1}{k}\right){\beta }_k^{\text{HS}}{d}_{k-1}^{\text{HS}}\right).\end{array}$ (9) The proposed method does not employ the hybrid conjugate gradient method [25] to generate the aperture shape. In contrast to the hybrid conjugate gradient method, the proposed method calculates and saves the corresponding $d_k^{\text{PRP}}$ and $d_k^{\text{HS}}$. The direction $d_k^{\text{PRP}-\text{HS}}$ is obtained by combining two conjugate gradients through the number of iterations used to determine the aperture shape. The aperture search direction $d_k^{\text{PRP}-\text{HS}}$ is not used to calculate $d_{k+1}^{\text{PRP}}$ and $d_{k+1}^{\text{HS}}$ in the next iteration. In the proposed method, only basic conjugate gradient classes are used to modulate the gradient. According to the Eq. (8), to generate an aperture shape in each iteration, PRP conjugate gradient descent direction $d_k^{\text{PRP}}$ and HS conjugate gradient descent direction $d_k^{\text{HS}}$ should be calculated according to the aperture gradient information and previous direction information. With an increase in k, the search direction of the aperture shape $d_k^{\text{PRP}-\text{HS}}$ gradually moves from the PRP-based conjugate gradient descent direction to the HS-based conjugate gradient descent direction. When the optimization search is close to the optimal value $\left(D_{js}^{*}\mathrm{,}{y}_{\kappa }^{*}\right)$, $g_{k-1}\approx g_k=0$, according to Eq. (7), the values of ${\beta }_k^{\text{PRP}}$ and ${\beta }_k^{\text{HS}}$ are approximately zero. The condition that must be satisfied for the optimal solution is as follows: ${\nabla }_{D_{js}\mathrm{,}{y}_{\kappa }}\mathcal{L}\left(D_{js}\mathrm{,}{y}_{\kappa }\mathrm{,}{\pi }_{js}\mathrm{,}{\rho }_{\kappa }\right)=\mathrm{0,}$ (10) where $\mathcal{L}\left(x\right)$ is the Lagrange function and ρκ is another Lagrange multiplier under the KKT conditions [5]. At $\left(D_{js}^{*}\mathrm{,}{y}_{\kappa }^{*}\right)$, subtract $\left(\frac{1}{k}{\beta }_{kjs}^{\text{PRP}}{d}_{\left(k-1\right)js}^{\text{PRP}}+\left(1-\frac{1}{k}\right){\beta }_{kjs}^{\text{HS}}{d}_{\left(k-1\right)js}^{\text{HS}}\right)$ from Eq. (10) without changing the optimal solution to the original optimization problem. The proposed method uses a conjugate gradient modulation to generate an aperture shape. Compared to generic column generation based on the negative gradient direction, the proposed method achieved faster convergence and shorter computation time for plan optimization.

Results from cases of head and neck cancer

Four optimization methods—generic column generation (labeled as "Original"), PRP conjugate gradient direction-based column generation (labeled as "PRP"), HS conjugate gradient direction-based column generation (labeled as "HS"), and the proposed method integrating PRP and HS conjugate gradient directions (labeled as "PRP-HS")—were employed to optimize case H1. The optimization results are shown in Fig. 3 and Table 3.

Fig. 3Full-screen View

(Color online) Optimization results of case H1. (a) is the DVH of the OARs; (b) is the DVH of the PTVs after optimization; (c), (d), (e), and (f) are the dose distribution of case H1 on the same surfaces after optimization of the four methods; (g) is the decrease in the objective function value in the iterative process; and (h) is the enlarged figure of (g).

The DV curves of the PTVs optimized by all the methods were mostly consistent(Fig. 3b). This conclusion can also be verified by the DV percentage, HI, and CN of the PTVs as shown in Table 3. It is apparent from the NTCP and gEUD of the OARs in Table 3 that PRP-HS can better protect the OARs while ensuring dose distribution to the PTVs. This conclusion is supported by the results shown in Figs. 3c-f. In particular, the parotid glands received the lowest dose after optimization of PRP-HS. Table 3, Fig. 3g, and h also demonstrate that the optimization time of PRP-HS was the shortest, and the decrease in the objective function value was the largest in the optimization iteration process.

In this study, the generic aperture shape generation is regarded as a search for the negative-gradient descent direction. To address the shortcomings of the negative-gradient descent direction, the search direction was constructed using the conjugate gradient direction. However, when the classical conjugate gradient direction was used to construct the search direction, a decrease in the objective function value was not ideal for the entire iterative process. The proposed method based on conjugate gradient modulation was used to determine the aperture shape. The proposed method improves the search direction of generic column generation. Subsequent cases were optimized using only these two methods to verify the improved efficacy of PRP-HS over generic column generation. For H2, the DVHs, dose distributions, and decreases in the objective function values obtained using these two methods are shown in Fig. 4. The performance details of the methods for case H2 are presented in Table 4.

Fig. 4Full-screen View

(Color online) Optimization results of case H2. (a) is the DVH of the OARs; (b) is the DVH of the PTVs; (c) and (d) are the dose distributions of case H2 on the same surfaces after optimization of the two methods; and (e) is the decrease in the objective function value in the iterative process.

Similar to case H1, case H2 shows that the dose distributions of the PTVs were mostly consistent (Fig. 4b and Table 4). Table 4 shows that, compared with the Original, PRP-HS can reduce the NTCP and gEUD of the OARs. However, the optimized mean dose received by the parotid glands and the maximum dose received by the spinal cord and brain stem did not satisfy the evaluation criteria listed in Table 1. The dose slices shown in Figs. 4c, d show that the organ distribution in case H2 was more compact than that in case H1, which makes it more difficult for the dose on the OARs, after optimization, to meet the evaluation criteria in Table 1. Fig. 4e shows that, compared with the Original, PRP-HS can significantly reduce the objective function value.

Results from cases of prostate cancer

Two optimization methods, the Original and PRP-HS methods, were used to optimize cases P1 and P2. Figure 5 depicts the optimization results, and Table 5 presents the evaluation index of the OARs and dose distribution to the PTV after optimization.

Fig. 5Full-screen View

(Color online) Optimization results of cases P1 and P2. (a) is the DVH for case P1; (b) is the DVH for case P2; (c) and (d) are the decreases in the objective function in the iterative process for cases P1 and P2 respectively; (e) and (f) are the dose distributions of case P1 on the same surfaces after optimization of the two methods; and (g) and (h) are the dose distributions of case P2 on the same surfaces after optimization of the two methods.

These two groups of comparative experiments verified the effectiveness of PRP-HS. In Figs. 5a, b, PRP-HS reduced the DV curves of the OARs in the high-dose region, while the DV curves of the PTV were mostly consistent. The PTV indicators in Table 5 confirm a similar dose distribution to the PTV, reducing the NTCP and gEUD of the OARs. Although the improvement in the plan quality was not as obvious as the optimization results of cases H1 and H2, the objective function value decreased significantly during the iterative process (Figs. 5c, d).

Supplementary experiments

The remaining cancer cases (labeled as “H3”, “H4”, “H5”, “P3”, “P4”, and “P5”) were used to verify the performance of the proposed method.

On the premise that the dose distribution to the PTVs is generally consistent, Tables 6 and 7 present the optimization information for the OARs. For the three head and neck cancer cases in Table 6, under the premise of meeting the dose constraints in Table 1, the proposed method reduces the dose received by the parotid glands. As shown in Table 7, the proposed method reduced the NTCP and gEUD of the OARs in the three cases of prostate cancer.

Discussion

A suitable gradient descent direction was selected for optimization to solve the pricing problem in generic column generation. This is considered the direction of negative gradient descent. Although it has a low computational complexity, the steepest descent method is prone to the sawtooth effect when searching near the minimum value, resulting in reduced search efficiency. Compared with Newton’s and quasi-Newton methods, the conjugate gradient method with its lower computational complexity is an ideal choice. Fig. 1 reveals that column generation methods based on conjugate gradient direction are superior to generic column generation in different degrees during the iteration. However, the convergence performance of column generation methods based on different conjugate gradient directions is not always suitable during the optimization process. A method based on conjugate gradient modulation was proposed to accelerate the descent speed of the objective function value without reducing the result quality. In Figs. 3g, h, the proposed method combined the advantages of PRP and HS, and decreased faster than the Original method. The optimization times listed in Table 3 show that PRP-HS was the fastest. Compared to column generation based on the negative gradient direction, in this optimization process, the column generation method based on conjugate gradient modulation made the decrease in the objective function value more evident (Figs. 4e, 5c, d).

When DVH is used for evaluation, the DV curves of the PTV optimized by all methods should be as similar as possible. Accordingly, the performance of the methods can be evaluated by observing the DVH of the OARs. For case H1, the DV curves of the PTVs in the DVH optimized using the four contrasting methods (Original, PRP, HS, and PRP-HS) were similar (Fig. 3b). The DVH of the OARs (Fig. 3a) showed that the DV curves obtained by PRP-HS were the lowest for the parotid glands. For the spinal cord and brain stem revealed that the DV curves optimized by PRP-HS were slightly worse than those obtained using the other three methods. However, the maximum doses received by the spinal cord and brain-stem are greater concern (Table 1). The maximum dose received by the spinal cord optimized using PRP and HS exceeded 50 Gy, whereas when optimized using PRP-HS, the result was minimal (Table 3). The maximum dose received by the brain stem optimized by the four methods satisfied $D_{\text{max}}\le 54$ Gy. For case H2, the gEUD of the spinal cord and brain stem optimized by PRP-HS were slightly worse than those obtained using the Original method (Table 4). By definition, the gEUD examines the entire dose distribution in the whole organ. As shown in Fig. 4a, the overall DV curves trend of the spinal cord and brain stem obtained with PRP-HS was worse than that of the Original method. According to Table 4, the maximum dose received by the brain stem optimized using both optimization methods satisfies $D_{\text{max}}\le 54$ Gy. The maximum dose received by the spinal cord, optimized by the Original dose, exceeded 50 Gy. Compared with the contrast method, PRP-HS reduced the NTCP and gEUD of the parotid glands to varying degrees and had a better ability to protect the OARs. These results illustrate that the proposed optimization method can better protect the OARs while ensuring dose distribution to the PTVs than the generic method. The optimization results for P1 and P2 (Fig. 5 and Table 5) also prove the proposed method performance.

For cases of head and neck cancer, only the evaluation index in Table 1 is presented in Table 6, whereas for cases of prostate cancer, only the NTCP and gEUD of the OARs are presented in Table 7. These results are sufficient to illustrate the improvement in the optimization for the proposed method compared with the generic method.

Table 8 presents the optimization results for the five cases of head and neck cancer, and the five cases of prostate cancer involved in this study were statistically analyzed. The proposed method significantly improved the dose distribution to the OARs (P < 0.05) other than the brain stem (P = 0.441), and the stability of the proposed method was proven.

Compared to generic column generation based on negative gradient direction, PRP-HS required less time to optimize the treatment plan. According to the analysis of the optimization results, compared with generic column generation, the proposed optimization method did not reduce and, in some cases, improved the plan quality. The proposed method did not reduce the number of generated apertures, which is a direction for future research and improvement. In the DAO, the aperture shape can be optimized by selecting the descending direction that maximizes improving the objective function. To generate a deliverable aperture shape, the descent direction was the direction after the conjugate gradient modulation in the pricing problem of this study. The proposed aperture shape optimization method does not involve line search or step length selection. In the proposed method, only basic conjugate gradient classes were used to modulate the gradient. In the future, some newly proposed conjugate gradient methods and radiation therapy techniques [38, 39] will be introduced into the DAO.