Li-Yuan Zhang · Zhi-Guo Gui · Peng-Cheng Zhang · Jie Yang

Abstract The gradient element of the aperture gradient map is utilized directly to generate the aperture shape without modulation.This process can be likened to choosing the direction of negative gradient descent for the generic aperture shape optimization.The negative gradient descent direction is more suitable under local conditions and has a slow convergence rate.To overcome these limitations, this study introduced conjugate gradients into aperture shape optimization based on gradient modulation.First, the aperture gradient map of the current beam was obtained for the proposed aperture shape optimization method, and the gradients of the aperture gradient map were modulated using conjugate gradients to form a modulated gradient map.The aperture shape was generated based on the modulated gradient map.The proposed optimization method does not change the optimal solution of the original optimization problem, but changes the iterative search direction when generating the aperture shape.The performance of the proposed method was verified using cases of head and neck cancer,and prostate cancer.The optimization results indicate that the proposed optimization method better protects the organs at risk and rapidly reduces the objective function value by ensuring a similar dose distribution to the planning target volume.Compared to the contrasting methods, the normal tissue complication probability obtained by the proposed optimization method decreased by up to 4.61%, and the optimization time of the proposed method decreased by 5.26% on average for ten cancer cases.The effectiveness and acceleration of the proposed method were verified through comparative experiments.According to the comparative experiments, the results indicate that the proposed optimization method is more suitable for clinical applications.It is feasible for the aperture shape optimization involving the proposed method.

Keywords Aperture shape · Column generation · Conjugate gradient · Gradient modulation · Direct aperture optimization

1 Introduction

Direct aperture optimization (DAO) [1] for intensity-modulated radiation therapy (IMRT) can be achieved by using approaches such as stochastic search, local gradient-based methods, and column generation.A stochastic search randomly moves the multileaf collimator (MLC) leaves to either side in small increments from their current position [2, 3].If the motion of the leaves improves the objective function, the current position is updated.Otherwise, the position change is accepted with a certain probability to avoid the local optima.The random nature of stochastic search makes it inefficient for aperture shape optimization.In the local gradient-based method, the positions of the MLC leaves are the optimization variables in the objective function of the optimization problem[4].Because this method uses the local gradient of positions to generate an aperture shape, it easily reaches the local optima.In practice, this method largely depends on an appropriate initial solution.An alternative is column generation, which notably differs from the other two methods [5, 6] that no initial aperture shape is required, and the gradient information is not local, as the network flow is constructed using the gradient map of the entire aperture to obtain a deliverable aperture shape.Unlike the two methods mentioned above, the optimization of column generation does not depend on the suitability of the initial solution.In this study, an improved method based on generic column generation is developed.

In one iteration of the generic method, the pricing problem is solved first, followed by the master problem.To solve the pricing problem, a network flow is constructed using the aperture gradient map.It is used to solve the shortest path problem and obtain a deliverable aperture shape.This process is equivalent to choosing a suitable descent direction for the objective function of the optimization problem.If the gradient element of the gradient map is not modulated, the process is equivalent to selecting the steepest descent direction for optimization.

The steepest descent method has the advantage of low computational cost and can converge from any initial point to a local minimum.However, this method typically exhibits a sawtooth effect in the region near the minimum value.Newton’s method achieves extremely fast convergence near the optimum, but is computationally expensive.Quasi-Newton methods avoid the explicit matrix required in conventional methods; however, the computation remains highly complex.The conjugate gradient method is an effective substitute because it has a comparable computational cost and converges faster than the steepest descent method.The computational complexity of the conjugate gradient method is less than that of the Newton’s and quasi-Newton methods.The conjugate gradient method is particularly suitable for solving large-scale optimization problems and is widely used in economics, engineering, physics, and other fields [7–9].In IMRT, the conjugate gradient has been used to optimize the weights of beams (apertures) [2, 10] and study the performance of column generation [11].

To overcome the drawbacks of aperture shape optimization based on the negative gradient descent direction, this study introduced conjugate gradients into aperture shape optimization.The search direction containing the conjugate gradient information was constructed using the original gradient of the gradient map to efficiently generate the aperture shape.

2 Methods

In this study, based on column generation, the aperture shape search direction containing the conjugate gradient components was constructed to overcome the problem of slow convergence in generating an aperture shape directly by using the gradient map without gradient modulation.

2.1 Dose calculations

During radiotherapy, the patient is irradiated with a predefined beam set denoted byB.Each beam in this study consists ofmrows andncolumns of beamlets, with each beamlet size being 1 cm × 1 cm.The set of all generated deliverable apertures is denoted byK, and the weight of the apertureκisyκ.The beamlets in setAκare delivered to the patient via apertureκ, and the treatment involvesSstructures, where each structures(s=1,…,S) comprisesvsvoxels.The element in the deposition matrix is the deposition coefficientWijs, which represents the dose received byj(j=1,…,vs) in structuresfrom the beamleti(i∈Aκ)allowed to pass through the apertureκat unit intensity.The doseDjscan be expressed as follows:

2.2 Column generation

In this study, the optimization problem is constructed as

When an aperture is generated, its cost should be calculated [5].The pricing problem is

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.

2.3 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.1 Nonlinear conjugate gradient method

Hestenes et al.[15] proposed a linear conjugate gradient method for solving linear equations, whereas Fletcher and Reeves [16] proposed a nonlinear conjugate gradient method for minimizing general functions.These two methods have been subsequently improved.In the (k+1)th iteration of the generic process of the nonlinear conjugate gradient method,the iteration format isxk+1=xk+αkdk.The conjugate gradient directiondkis updated as follows:

wherexdenotes the independent variable of the objective function,ddenotes the search direction (i.e., conjugate gradient direction), andgdenotes the first derivative of the objective function.The parameterβkis important for nonlinear conjugate gradient methods, because it determines the type of method used.Since 1952, a series of representative conjugate gradient methods have been proposed, such as the Fletcher–Reeves (FR) [16], Polak–Ribère–Polyak (PRP)[17, 18], Hestens–Stiefel (HS) [15], Dai–Yuan (DY) [19],conjugate descent (CD) [20], and Liu–Storey (LS) [21] conjugate gradient methods.The common definitions ofβkare as follows.

where ‖?‖ represents the Euclidean norm andyk-1=gk-gk-1.

In recent years, many scholars have improved the conjugate coefficientβfor the conjugate gradient method [22] to achieve good convergence and numerical results [23, 24].This study attempted the conjugate gradient modulation for gradient elements of the aperture gradient map.Only the basic conjugate gradient direction was used at this stage.

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).

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.

3 Experimental setup and evaluation criteria

In this experiment, the classical pencil beam method [26] in an open-source computational environment for radiotherapy research (CERR) [27] was used to calculate the dose deposition matrixW.All methods involved in the experiment were implemented in Visual C++ (version Visual Studio 2012)on a computer with an Intel?CoreTMi9-10900X central processing unit at 3.70 GHz, running on Windows 10 with 64 bits.All cancer cases involved in this study were obtained from Shanxi Provincial Cancer Hospital.The simulation experiments involving the relevant cancer cases in this study were conducted under a protocol approved by the Ethics Committee of the North University of China.

The physician defined the organs involved in the head and neck cancer cases, as illustrated in Fig.2a, based on image data.Three planning target volumes (PTVs), PTV 70 Gy,PTV 63 Gy, and PTV 56 Gy, were obtained by an outward expansion of 5 mm in each direction of the three clinical target volumes (CTVs).The ipsilateral parotid gland (IL-PG),contralateral parotid gland (CL-PG), spinal cord, and brain stem were considered as the organs at risk (OARs) in the optimization [28].Those OARs were obtained by expanding 5 mm outwards from the outline of each organ.Nine equally spaced 6 MeV co-irradiated photon fields were set on the CERR to simulate the head and neck cancer case irradiation.In the optimization process, we constrained the dose distribution to the parotid glands using the dose–volume(DV) criterion [29] sub-objective function and to the spinal cord and brain stem using the maximum dose criterion to penalize doses beyond the upper limit, as listed in Table 1.The dose distributions to the PTVs were constrained using the mean and minimum dose sub-objective functions, and the dose distribution to the remaining normal tissues was constrained using the maximum dose sub-objective function.The optimization iterations were capped at 100 iterations.

Fig.2 (Color online) Distribution of organs involved in the cancer cases: a for the head and neck cancer case and b for the prostate cancer case

For the prostate cancer cases illustrated in Fig.2b, one PTV was set, and the OARs were the bladder and rectum.The OARs and PTV were delineated by the physician using the image data.Manually delineated bladder and rectal contours were expanded outwards by 5 mm to obtain the OARs.The CTV was expanded backward by 5 mm and outwards by 10 mm in the remaining directions to obtain the PTV.Five radioactive sources with gantry angles of 36?, 100?,180?, 260?, and 324?were used for simulated irradiation.The dose distributions for the bladder and rectum were constrained using the DV criterion sub-objective function, and the DV constraint conditions are listed in Table 2.The mean and minimum dose criteria constrained the dose distribution to the PTV.The maximum dose sub-objective function restrained the remaining normal tissue dose distribution.The iterations of the optimization process were up to 60.

4 Results and discussion

Comparative experiments were conducted to compare the performance of different methods for head and neck cancer cases (labeled as “H1”, “H2”, “H3”, “H4”, and “H5”) and prostate cancer cases (labeled as “P1”, “P2”, “P3”, “P4”, and“P5”).The total objective function used in the experiments was the sum of multiple sub-objective functions multiplied by the corresponding penalty factors [30].In the same set of experiments, different contrast methods used the same objective function and penalty factors of the sub-objective functions.In Sects.4.1 and 4.2, the optimized results for cases H1, H2, P1, and P2 are presented.To evaluate the optimized results for cases H1, H2, P1, and P2, the DV histogram (DVH) was analyzed using the clinical guidance standard (Tables 1, 2) developed by Marks et al.[31].The generalized equivalent uniform dose (gEUD) and normal tissue complication probability (NTCP) of head and neck[32, 33] and prostate cancer cases [34, 35] were calculated to evaluate the protective effect of each method on the OARs.Relatively lower gEUD and NTCP values indicated better protection of the OARs.The conformity number (CN) [36]and homogeneity index (HI) [37] of the PTV were calculated.When CN and HI were close to 1, the dose distribution to the PTV was more conformal and uniform.The running time, number of apertures, and trend of the objective function during optimization were also used to investigate the performance of the experimental methods.The optimized results for the remaining six cancer cases are concisely presented.

Table 1 DV constraint conditions for the head and neck cancer case

Table 2 DV constraint conditions for the prostate cancer case

4.1 Results from cases of head and neck cancer

Four optimization methods–generic column generation(labeled as “Original”), PRP conjugate gradient directionbased 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.

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 Fig 3c–f.In particular, the parotid glands received the lowest dose after optimization of PRP-HS.Table 3, Fig.3g, 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 case 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.

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 Fig.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.Figure 4e shows that, compared with the Original, PRP-HS can significantly reduce the objective function value.

4.2 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.3 (Color online) Optimization results of case H1.a DVH of the OARs; b DVH of the PTVs after optimization; c–f dose distributions of case H1 on the same surfaces after optimization of the four methods; g decrease in the objective function value in the iterative process; and h enlarged figure of g

Table 3 Results obtained by the four optimization methods for case H1

These two groups of comparative experiments verified the effectiveness of PRP-HS.In Fig.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 (Fig.5c, d).

4.3 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.

4.4 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.Figure 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 Fig.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).

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

Table 4 Results obtained by the two optimization methods for case H2

Table 5 Results obtained by the two optimization methods for cases P1 and P2

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 PRPHS) 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, it was revealed that the DV curves optimized by PRPHS 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 satisfiedDmax≤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 curve 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 satisfiesDmax≤54 Gy.The maximum dose received by the spinal cord, optimized by the Original, 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 proved.

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.

5 Conclusion

This study proposed a method based on conjugate gradient modulation for aperture shape optimization by modulating the gradients, because the negative gradient direction search exhibited a sawtooth effect when approaching the local optimal point.The conjugate gradient method with low computational cost and fast convergence was introduced into the aperture generation.The performance of PRP-HS was verified in head and neck cancer cases, and prostate cancer cases.Based on comparative experimental results,the proposed optimization method can accelerate the solution process and improve the quality of the treatment plan.The optimization time of the proposed method decreased by up to 9.41% for head and neck cancer cases, and 9.71%for prostate cancer cases.This improvement does not come at the expense of the quality of the results.While ensuring dose distribution to the PTV, PRP-HS reduced NTCP by up to 4.61% compared with generic column generation.According to the experimental results, the proposed aperture shape optimization method can be applied to radiotherapy plan optimization for different cancer cases and can be effi-ciently used in clinical settings.

Table 6 Results for cases H3,H4, and H5

Table 7 Results for cases P3,P4, and P5

Table 8 Statistical analysis of the optimization results

Author contributionsAll authors contributed to the study conception and design.Material preparation, data collection, and analysis were performed by Li-Yuan Zhang, Zhi-Guo Gui, Peng-Cheng Zhang,and Jie Yang.The first draft of the manuscript was written by Li-Yuan Zhang, and all authors commented on previous versions of the manuscript.All authors read and approved the final manuscript.

Data availabilityThe data that support the findings of this study are openly available in Science Data Bank at https:// www.doi.org/ 10.57760/ scien cedb.j00186.00171 and https:// cstr.cn/ 31253.11.scien cedb.j00186.00171.

Declarations

Conflict of interestThe authors declare that they have no competing interests.