An Inverse Dose Optimization Algorithm for Three-Dimensional Brachytherapy

Purpose: To investigate an implementation method and the results of an inverse dose optimization algorithm, Gradient Based Planning Optimization (GBPO), for three-dimensional brachytherapy.

Methods: The GBPO used a quadratic objective function, and a dwell time modulation item was added to the objective function to restrict the dwell time variance. We retrospectively studied 4 cervical cancer patients using different applicators and 15 cervical cancer patients using the Fletcher applicator. We assessed the plan quality of GBPO by isodose lines for the patients using different applicators. For the 15 patients using the Fletcher applicator, we utilized dose-volume histogram (DVH) parameters of HR-CTV (D_100%, V_150%) and organs at risk (OARs) (D_0.1cc, D_1cc, D_2cc) to evaluate the difference between the GBPO plans and the IPSA (Inverse Planning Simulated Annealing) plans, as well as the GBPO plans and the Graphic plans.

Results: For the 4 patients using different applicators, the dose distributions are conformable. For the 15 patients using the Fletcher applicator, when the dwell time modulation factor (DTMF) is less than 20, the dwell time deviation reduces quickly; however, after the DTMF increased to 100, the dwell time deviation has no remarkable change. The difference in dosimetric parameters between the GBPO plans and the IPSA plans is not statistically significant (P>0.05). The GBPO plans have a higher D_100% (3.57 ± 0.36, 3.38 ± 0.34; P<0.01) and a lower V_150% (55.73 ± 4.06, 57.75 ± 3.79; P<0.01) than those of the Graphic plans. The differences in other DVH parameters are negligible between the GBPO plans and the Graphic plans.

Conclusions: The GBPO plans have a comparable quality as the IPSA plans and the Graphic plans for the studied cervical cancer cases. The GBPO algorithm could be integrated into a three-dimensional brachytherapy treatment planning system after studying more sites.

Introduction

Compared with external beam radiation therapy, brachytherapy has the characteristics of high dose near the source and rapid dose drop-off away from the source. In addition, because the applicator is implanted in the tumor region, brachytherapy reduces the dosimetric uncertainties caused by anatomical change and setup error. These advantages ensure the irreplaceable role of brachytherapy in radiotherapy (1).

At the present stage, image-guided three-dimensional (3D) brachytherapy is the mainstream method for brachytherapy. Dose optimization is a crucially important component of 3D brachytherapy treatment planning systems (TPSs). In general, dose optimization methods of 3D brachytherapy can be divided into forward optimization and inverse optimization. In a forward optimization process, a planner manually enters the dwelling weight/time or drags isodose lines based on the planner’s clinical experience to achieve a desirable dose distribution. The method of dragging isodose lines is called graphic optimization. In an inverse optimization process, a planner inputs the objectives and penalty weights of targets and organs at risk (OARs) based on the prescription dose and patient’s anatomy. Through a trial-and-error process, a satisfactory dose distribution can be generated by the inverse dose optimization system. Inverse optimization algorithms of brachytherapy, such as IPSA (Inverse Planning Simulated Annealing) and HIPO (Hybrid Inverse Planning Optimization), have been reported in literatures and implemented in 3D TPSs (2–5).

The inverse optimization algorithm of brachytherapy usually produces a plan with a large dwell time variation (6), which should be addressed for the following reasons: First, a location with a large dwell time is suspect to have a high dose. A high dose region should be avoided unless a tumor volume requires an inhomogeneous dose distribution. Second, the larger the dwell time variation, the greater the inhomogeneous dose distribution. An inhomogeneous dose distribution is more likely to be affected by source position uncertainties. Both IPSA and HIPO provide parameters that restrict dwell time variance: the Dwell Time Deviation Constraint (DTDC) and Dwell Time Gradient Restriction (DTGR) for IPSA and HIPO, respectively (4, 5). By adjusting these parameters, it is possible to obtain a favorable clinical plan for which the variation of the dwell time is considered acceptable.

The purpose of this paper is to investigate an in-house inverse brachytherapy optimization algorithm referred to as Gradient Based Planning Optimization or GBPO and a new method to restrict dwell time variance. We retrospectively studied a total of 4 cervical cancer patients using different applicators and 15 cervical cancer patients using the Fletcher applicator to evaluate the GBPO algorithm.

Materials and Methods

Dose Calculation

The dose calculation algorithm in this study was based on the AAPM TG-43 recommendation (7, 8). Since the implementation detail has been reported in the reference (9), only a brief introduction is included here. We calculated the dose of the i-th voxel, D_i, through the formula given in Equation 1:

$\begin{array}{ll}{D}_i={\displaystyle {\sum }_{m=1}^{N_M}{\displaystyle {\sum }_{N=1}^{N_N}{d}_{m,n}{t}_{m,n}}}& (1)\end{array}$

where N_M is the total channel number, N_N is the total dwell position number in the m-th channel, d_m,n and t_m,n are the dose rate contribution and the dwell weight, respectively, from the n-th dwell position in the m-th channel.

Inverse Optimization

The GBPO optimization algorithm was implemented using the LBFGS (Limited memory Broyden Fletcher Goldberg Shanno) code, which is an optimization engine based on the gradient descent method (10, 11). The GBPO used a quadratic objective function, and we calculated the objective value F through the formula given in Equation 2:

$\begin{array}{ll}\begin{array}{l}F\left(t_{m,n}\right)={{\displaystyle {\sum }_{i\in TAR}{p}_{TAR}·H_{TAR}\left(D_i-D_{0,TAR}\right)·\left(D_i-D_{0,TAR}\right)}}^2\hfill \\ +{\displaystyle {\sum }_{i\in OARs}{p}_{OARs}·H_{OARs}(D_i-}{D}_{0,OARs})·{\left(D_i-D_{0,OARs}\right)}^2\hfill \\ +p_{SOU}·{\displaystyle {\sum }_{m=1}^{N_M}{\displaystyle {\sum }_{n=1}^{N_N}\frac{1}{N_N}{\left(t_{m,n}-t_{m,\min }\right)}^2}}\hfill \end{array}& (2)\end{array}$

where p_TAR is the penalty weight of the target; H (D_i – D₀) is a Heaviside function (12), and for a target it equals 0 if D_i > D₀ but 1 if D_i ≤ D₀; the value reverses for an OAR. D_i is the dose of the i-th voxel; D₀ is the objective dose; p_TAR is the penalty weight of the OARs; p_SOU is the dwell time modulation factor (DTMF); and t_m,min is the smallest dwell time in the m-th channel.

The GBPO considered multiple targets and OARs, and each region of interest (ROI) had an objective dose. The last item in Equation (2) was provided to modulate the dwell time variance to meet the clinical needs.

Test Cases

We divided the clinical test of the GBPO algorithm into two parts: the first part tested the optimization results of different applicators, which include a double ovoid applicator, a tandem-ring applicator, a multi-channel applicator, and a tandem-needles applicator. In the second part, we retrospectively studied 15 cervical cancer patients using the Fletcher applicator, and the average HR-CTV volume was 52.65 cm³ (minimum 36.03 cm³, maximum 80.45 cm³).

Treatment Planning

The delineation of target and OARs was performed on an Oncentra V4.3 (Elekta AB, Stockholm, Sweden) TPS. The target was HR-CTV, and the OARs included the bladder, rectum, and sigmoid. The dose prescription was 6 Gy.

For each patient, we compared the following three plans: the IPSA plan, the Graphic plan, and the GBPO plan. For all plans, the source step size was 0.25 cm, and the dose calculation grid resolution was 0.1 cm x 0.1 cm x 0.1 cm. Since the optimization results were affected by the dwell point number and dwell position, the three plans used the same dwell point number and dwell position.

The IPSA plan automatically determined dwell positions based on the reference target. The DTDC value affects the optimization result (13), and we set it to 0.4 for all IPSA plans in this study.

We changed the dwell time of each dwell position to 1 before the Graphic optimization, and then a physicist manually dragged the isodose line to achieve a desirable dose distribution. The quality of the Graphic plan heavily depends on the clinical experience of the physicist. In order to improve the quality of the Graphic plan, the planning was performed by an experienced physicist who has worked in the brachytherapy department for more than 5 years.

For the GBPO plan, we set the initial dwell time of each dwell position to 1. The minimum value of the dwell time in the GBPO iteration process was set to 0.000001, which ensures the non-negativity of the dwell time during the optimization process. All GBPO plans were iterated 100 times. Table 1 gives the optimization objectives for studying the relationship between the DTMF and the dwell time standard deviation (DTSD), as well as the initial objectives of the GBPO plans used for the comparison with the IPSA plans. In the dosimetric comparison process, if the optimization result of a GBPO plan was not satisfactory, we adjusted the initial objectives until obtaining a satisfying result. We set the DTMF to 10 for all GBPO plans, based on the results given in Figure 2.

TABLE 1

Table 1 The optimization objectives for GBPO plans (Dose unit: Gy).

Plan Evaluation

We assessed the plan quality of GBPO by isodose lines for the 4 patients using different applicators. For the 15 patients using the Fletcher applicator, dose volume histogram (DVH) parameters were used to evaluate the dosimetric difference between the GBPO plans and the IPSA plans, as well as the GBPO plans and the Graphic plans. We defined D_x% as the dose expressed in Gy that received by x% of the total volume, V_y% as the volume expressed in percentage that received y% of the prescribed dose, and D_zcc as the dose expressed in Gy that received by z cm³ volume. The DVH parameters for HR-CTV were D_100% and V_150%, and for the OARs, they were D_0.1cc、D_1ccand D_2cc (14). All plans were normalized to HR-CTV D_90% =6Gy. To evaluate the dosimetric parameters mentioned above, the SPSS 19.0 software (IBM, Armonk, NY) was used for statistical analyses. We conducted paired, a two-sided Wilcoxon signed-rank test to compare the dose distributions between the GBPO plans and the IPSA plans, as well as the GBPO plans and the Graphic plans. A P value less than 0.05 was considered statistically significant.

Results

Isodose Lines

Figure 1 shows the isodose lines optimized by the GBPO algorithm for the 4 patients using different applicators. The GBPO algorithm can generate conformable dose distributions for different applicators.

FIGURE 1

Figure 1 Isodose line of patients using different applicators. (A) Double ovoid applicator; (B) Tandem-ring applicator; (C) Multi-channel applicator; (D) Tandem-needles. The organs are HR-CTV (red), bladder (blue), rectum (brown), and sigmoid (green).

Dwell Time Modulation Factor

Figure 2 illustrates the DTSD of the 15 patients using the Fletcher applicator optimized by the GBPO algorithm. When the DTMF is less than 20, the DTSD decreases quickly, but the DTSD has no remarkable change after the DTMF increased above 100. Therefore, in the planning optimization of cervical cancer, a DTMF value greater than 100 is not recommended when using the GBPO algorithm.

FIGURE 2

Figure 2 The dwell time standard deviation (DTSD) as a function of dwell time modulation factor (DTMF) for patients using the Fletcher applicator.

Dosimetric Comparison

Table 2 compares the DVH parameters of the target and OARs for the 15 patients. The difference in dosimetric parameters between the GBPO plans and the IPSA plans is not statistically significant (P > 0.05). The GBPO attains a similar plan quality as the IPSA. The GBPO plans has a higher D_100% (3.57 ± 0.36, 3.38 ± 0.34; P<0.01) and a lower V_150% (55.73 ± 4.06, 57.75 ± 3.79; P<0.01) than that of the Graphic plans. The differences in other DVH parameters are negligible between the GBPO plans and the Graphic plans.

TABLE 2

Table 2 A dosimetric comparison for target and organs at risk.

Discussion

With the aid of imaging techniques such as CT and MRI, we can obtain an applicator position and ROIs three dimensionally. Knowing the applicator position and the source step size, we can also determine the dwell positions. A variety of algorithms have been developed to optimize the dwell time to achieve a desirable dose distribution (2–4, 15–18)_ENREF_11. In this study, we implemented a new inverse optimization algorithm, GBPO, to optimize the 3D brachytherapy dose based on patient anatomy and prescription dose. The patient data show that this algorithm achieves similar optimization results as compared with a commercial algorithm.

Uncertainties affect dose accuracy in high dose-rate brachytherapy (19). Regional hotspots should be avoided. Several studies have suggested using the dwell time modulation factor to address the issue of large dwell time variation (20, 21). In IPSA, the DTDC is a user-entered parameter that constrains the upper limit of a single dwell time relative to the average dwell time (4). The DTDC changes in the [0-1.0] range by a step of 0.1. When the DTDC is 0, it means that the optimization has no dwell time constraint, and the dwell time is the most uniform when it equals 1.0. The DTDC effectively reduces large dwell times, and neglects dwell times below the average value. The dwell time modulation factor of HIPO is DTGR, which avoids a large dwell time change between adjacent dwell locations, and eliminates the existence of large dwell times that may cause hotspots. Similar to the DTDC, the DTGR varies by a step size of 0.1 in the range of (0–1.0) (5). Increasing the DTGR value forces the optimizer to avoid situations where the dwell time is very long or very short. Since the DTGR considers the change of the adjacent dwell time, in places where there is no need to dwell, there may also be short dwell times if using DTGR.

The dwell time modulation principle of this algorithm is different from HIPO and IPSA. First, GBPO used the minimum dwell time in the objective function instead of the average value. The reason for this is that some of the dwell positions may not be suitable for dwell due to the OAR’s constraint, in which case the minimum dwell time is retained in GBPO. The purpose of adding a dwell time modulation item to the objective function is to make the larger dwell time shorter. Since the minimum dwell time is preserved, the DTSD will not be zero, even if the DTMF increases. Second, the GBPO does not normalize the maximum DTMF to 1. This study used the site of radical cervical cancer for testing, and did not consider other cancer sites. The normalized DTMF for one site may not be suitable for another site, so there are limitations in its application for other sites. In addition, using the same normalization method as IPSA and HIPO will make the modulation space smaller. There are only 11 values after normalization. Without normalization, we can change the DTMF value according to different clinical requirements.

There is currently a high incidence of cervical cancer (22). External beam radiation therapy combined with brachytherapy is the standard radiotherapy mode for cervical cancer (23). In our brachytherapy center, more than 90% patients have cervical cancer, which is why we selected cervical cancer patients as test cases. The Fletcher applicator is one of the most commonly used applicators for radical cervical cancer cases, it has 3 channels, and the optimization freedom is limited. However, compared with the IPSA plan, the GBPO obtained more favorable results, which gives us confidence this algorithm could be extended to other applicators and tumor sites.

For patients with radical cervical cancer, when the DTMF exceeds 100, the change in DTSD is not remarkable (Figure 1). Therefore, it is recommended to select the DTMF value within (0–100) for radial cervical cancers. It should be noted that Figure 2 is based on the optimization parameters listed in Table 1. The relation of the DTSD and the DTMF may vary if the optimization parameter changes. Different sites may have a different DTMF- DTSD curve. Therefore, the DTMF- DTSD curve of other sites should be studied before using DTMF to determine an appropriate value suitable for other cancer sites. Testing the applicability of this algorithm to other cancers is a topic of our future work.

Dose optimization is a trial-and-error process, for a reverse optimization algorithm, the calculation speed is a factor that needs to be considered. The GBPO calculation is performed on a single central processor unit now, so the time required for GBPO is longer than that of the IPSA for the same test case. The GBPO running time for each case is about 2 to 5 min, depending on the quantity of dwell positions and dose points. In order to reduce the time spent on dose optimization, the GBPO needs parallel computing by MPI (Message Passing Interface) or CUDA (Compute Unified Device Architecture) technique. Parallel computing is the work we are currently doing.

Conclusions

In this paper, we investigated a new inverse optimization algorithm, GBPO, for 3D brachytherapy, including a new dwell time modulation method. For a commonly used applicator in cervical cancer, this algorithm achieved similar results as compared with the IPSA optimization method. The GBPO algorithm could be integrated into a 3D brachytherapy TPS after more cancer sites are studied.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material; further inquiries can be directed to the corresponding author.

Ethics Statement

The studies involving human participants were reviewed and approved by Sichuan Cancer Hospital Ethics Committee. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author Contributions

XW, PW, JLi: Study design. XW, BT, SK, LL, JLa: Data collection and data analysis. XW, PW, LO, JLi: Manuscript preparation, manuscript editing. QH, ZW, CG: Critical discussion of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by National Key Research and Development Project (No. 2016YFC0105103, No. 2017YFC0113100) and Chengdu Science and Technology Project (2019-YF09-00095-SN).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Ke Li from the University of Wisconsin for editing assistance.

References

1. Tanderup K, Menard C, Polgar C, Lindegaard J, Kirisits C, Pötter R. Advancements in brachytherapy. Adv Drug Deliver Rev (2017) 109:15–25. doi: 10.1016/j.addr.2016.09.002

CrossRef Full Text | Google Scholar

2. Lessard E, Pouliot J. Inverse planning anatomy-based dose optimization for HDR-brachytherapy of the prostate using fast simulated annealing algorithm and dedicated objective function. Med Phys (2001) 28(5):773–9. doi: 10.1118/1.1368127

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Lahanas M, Baltas DN. A hybrid evolutionary algorithm for multi-objective anatomy-based dose optimization in high-dose-rate brachytherapy. Phys Med Biol (2003) 48(3):399. doi: 10.1088/0031-9155/48/3/309

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Cunha A, Siauw T, Hsu IC, Pouliot J. A method for restricting intracatheter dwell time variance inhigh-dose-rate brachytherapy plan optimization. Brachytherapy (2016) 15(2):246–51. doi: 10.1016/j.brachy.2015.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Gorissen BL, Dick DH, Hoffmann AL. Mixed integer programming improves comprehensibility and plan quality in inverse optimization of prostate HDR brachytherapy. Phys Med Biol (2013) 58(4):1041–57. doi: 10.1088/0031-9155/58/4/1041

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Chajon E, Dumas I, Touleimat M, Magné N, Coulot J, Verstraet R, et al. Inverse planning approach for 3-D MRI-based pulse-dose rate intracavitary brachytherapy in cervix cancer. Int J Radiat Oncol Biol Phys (2007) 69(3):955–61. doi: 10.1016/j.ijrobp.2007.07.2321

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Nirogi R, Kandikere V, Bhyrapuneni G, Benade V, Mekala V, Jayarajan P, et al. Dosimetry of interstitial brachytherapy sources: Recommendation of the AAPM Radiation Therapy Committee Task Group No. 43. Med Phys (1995) 22(8):209–34. doi: 10.1118/1.597556

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Rivard MJ, Coursey BM, Dewerd LA, Hanson WF, Huq MS, Ibbott GS, et al. Update of AAPM Task Group No. 43 Report: A revised AAPM protocol for brachytherapy dose calculations. Med Phys (2004) 31(3):633–74. doi: 10.1118/1.1646040

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Wang X, Wang P, Li C, Zhangwen W, Chengjun G, Li J, et al. An automated dose verification software for brachytherapy. J Contemp Brachytherapy (2018) 10(5):478–82. doi: 10.5114/jcb.2018.79396

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Liu DC, Nocedal J. On the limited memory BFGS method for large scale optimization. Math Programming (1989) 45(1-3):503–28. doi: 10.1007/BF01589116

CrossRef Full Text | Google Scholar

11. Byrd RH, Lu P, Nocedal J. Zhu C. A Limited Memory Algorithm for Bound Constrained Optimization. Siam J Sci Comput (1995) 16(5):1190–208. doi: 10.1137/0916069

CrossRef Full Text | Google Scholar

12. Bracewell R. Heaviside’s Unit Step Function. The Fourier Transform and Its Applications, 3rd ed. New York: McGraw-Hill (2000). p. 61–5.

Google Scholar

13. Smith RL, Panettieri V, Lancaster C, Mason N, Franich RD, Millar JL. The influence of the dwell time deviation constraint (DTDC) parameter on dosimetry with IPSA optimisation for HDR prostate brachytherapy. Australas Phys Eng Sci Med (2015) 38(1):55–61. doi: 10.1007/s13246-014-0317-2

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Potter R, Haie-Meder C, Van Limbergen E, Barillot I, De Brabandere M, Dimopoulos J, et al. Recommendations from gynaecological (GYN) GEC ESTRO working group (II): concepts and terms in 3D image-based treatment planning in cervix cancer brachytherapy-3D dose volume parameters and aspects of 3D image-based anatomy, radiation physics, radiobiology. Radiother Oncol (2006) 78(1):67–77. doi: 10.1016/j.radonc.2005.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Cui S, Despres P, Beaulieu L. A multi-criteria optimization approach for HDR prostate brachytherapy: I. Pareto surface approximation. Phys Med Biol (2018) 63(20):205004. doi: 10.1088/1361-6560/aae24c

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Cui S, Despres P, Beaulieu L. A multi-criteria optimization approach for HDR prostate brachytherapy: II. Benchmark against clinical plans. Phys Med Biol (2018) 63(20):205005. doi: 10.1088/1361-6560/aae24f

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Nicolae A, Morton G, Chung H, Loblaw A, Jain S, Mitchell D, et al. Evaluation of a Machine-Learning Algorithm for Treatment Planning in Prostate Low-Dose-Rate Brachytherapy. Int J Radiat Oncol Biol Phys (2017) 97(4):822–9. doi: 10.1016/j.ijrobp.2016.11.036

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Zhou Y, Klages P, Tan J, Chi Y, Stojadinovic S, Yang M, et al. Automated high-dose rate brachytherapy treatment planning for a single-channel vaginal cylinder applicator. Phys Med Biol (2017) 62(11):4361–74. doi: 10.1088/1361-6560/aa637e

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Kirisits C, Rivard M, Baltas D. Review of clinical brachytherapy uncertainties: analysis guidelines of GEC-ESTRO and the AAPM. Radiother Oncol (2014) 110(1):199–212. doi: 10.1016/j.radonc.2013.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Asa H, Torbjorn L, Asa CT. Impact of using linear optimization models in dose planning for HDR brachytherapy. Med Phys (2012) 39(2):1021–8. doi: 10.1118/1.3676179

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Marleen B, Gorissen BL, Dick DH, Hoffmann AL. Dwell time modulation restrictions do not necessarily improve treatment plan quality for prostate HDR brachytherapy. Phys Med Biol (2015) 60(2):537–48. doi: 10.1088/0031-9155/60/2/537

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: A Cancer J Clinicians (2018) 0(0):1–31. doi: 10.3322/caac.21492

CrossRef Full Text | Google Scholar

23. Wang X, Li J, Wang P, Yuan K, Yin G, Wan B. Image guided radiation therapy boost in combination with high-dose-rate intracavitary brachytherapy for the treatment of cervical cancer. J Contemp Brachytherapy (2016) 8(2):122–7. doi: 10.5114/jcb.2016.59282

PubMed Abstract | CrossRef Full Text | Google Scholar