Abstract
Charged-particle therapy is a rapidly growing precision radiotherapy technique that treats tumors with ion beams. Because ion-beam delivery systems have multiple degrees of freedom (including the beam trajectories, energies and fluences), it can be extremely difficult to find a treatment plan that accurately matches the dose prescribed to the tumor while sparing nearby healthy structures. This inverse problem is called inverse treatment planning (ITP). Many ITP approaches have been proposed for the simpler case of X-ray therapy, but the work dedicated to charged-particle therapy is usually limited to optimizing the beam fluences given the trajectories and energies. To fill this gap, we consider the problem of simultaneously optimizing the beam trajectories, energies, and fluences, which we call full ITP. The solutions are the global minima of an objective function defined on a very large search space and having deep local basins of attraction; because of this difficulty, full ITP has not been studied (except in preliminary work of ours). We provide a proof of concept for full ITP by showing that it can be solved efficiently using simulated annealing (SA). The core of our work is the incremental design of a state exploration mechanism that substantially speeds up SA without altering its global convergence properties. We also propose an original approach to tuning the cooling schedule, a task critical to the performance of SA. Experiments with different irradiation configurations and increasingly sophisticated SA algorithms demonstrate the benefits and potential of the proposed methodology, opening new horizons to charged-particle therapy.
Similar content being viewed by others
Notes
The local minima and their basins depend on how the search space is explored. By claiming that the objective has deep local basins, we actually mean that a landscape defined by the objective and a tractable exploration mechanism will likely have many deep (possibly poor) local minima. This fact will be confirmed by our experiments.
In fact, the full ITP problem is an instance of the so-called minimum penalty treatment problem, which is APX-complete [60].
The number of active orientations could be controlled by adding a group-sparsity penalty to the objective, as in [33], but this topic is beyond the scope of this paper.
Note that \(\Theta _1\) has the form (3.3) with weights \(\theta (\varvec{\omega },\varvec{\omega }') = |\mathcal {E}(\tau )|^{-1}\).
Although the possibility cannot be excluded, the argument of the minimum in (4.3) is unlikely to contain more than one element; so we assume it is a singleton for simplicity.
Therefore \(\{\mathcal {S}_2(\varvec{\omega })\}_{\varvec{\omega }\in \Omega } \) is not a neighborhood system.
Again, for simplicity, we assume that the argument of the minimum is a singleton.
The sequence \((p_m)_{m\in [1..M]}\) is also piecewise constant, with a number of plateaus less than or equal to the number of temperature stages.
This approximation is justified for charged particles, with a greater accuracy for heavy ions such as carbon than for light ions such as protons.
Throughout the paper, flops is the plural of flop (an elementary floating point operation) and is not to be confused with “flops per second”.
We recall that \(D(v,\varvec{b}'_\ell )\) is the product of the fluence of \(\varvec{b}'_\ell \) and the fluence-normalized dose \(D(v,(\tau ,e,1))\).
Most of the time, \(\mathcal {V}(\varvec{b}_\ell ) \cap \mathcal {V}(\varvec{b}'_\ell ) = \emptyset \) or \(|\mathcal {V}(\varvec{b}_\ell ) \cap \mathcal {V}(\varvec{b}'_\ell ) | \ll |\mathcal {V}(\varvec{b}_\ell )| + |\mathcal {V}(\varvec{b}'_\ell )|\) (unless the center axes of \(\varvec{b}_\ell \) and \(\varvec{b}'_\ell \) are parallel and close to each other), and so \(N_1(\varvec{\omega },\varvec{\omega }') \approx 4|\mathcal {V}(\varvec{b}_\ell )| + 5|\mathcal {V}(\varvec{b}'_\ell )|\).
This observation does not apply to SA with logarithmic cooling.
The point \((\log T_l, \log {\widehat{\zeta }}_{l})\) such that \({\widehat{\zeta }}_{l-1} \geqslant \frac{1}{2} \) and \({\widehat{\zeta }}_{l} < \frac{1}{2} \) must then be discarded from the regression data.
The total number of distinct beam trajectories in such a plan is therefore at most 4800, which remains much smaller than the number of possible trajectories \(|\mathcal {T}|\).
The subscript “opt” in \({\mathbb {A}}_{3, {\mathrm {opt}}}\) stands for “optimized” in reference to the incremental design of the proposed exploration strategy, and “p.c.” stands for “piecewise constant”
Because \(\mathcal {S}_{3,p}(\varvec{\omega }) \subset \mathcal {S}_2(\varvec{\omega }) \subset \mathcal {S}_1(\varvec{\omega })\) for all \(p\in {\mathbb {N}}^*\) and \(\varvec{\omega }\in \Omega \).
Plotting \(U(\varvec{\omega }_M)\) versus M is of little interest because the different communication mechanisms have different computational costs (see Sect. 5) .
When \({\mathscr {D}}\) is not stored in memory, the CPU time is multiplied by about \(7 \times R^{-1} |\mathcal {I}|\) for \({\mathbb {A}}_{1 , \log }\) and \({\mathbb {A}}_{1 , \exp }\), \(4.8 \times R^{-1} |\mathcal {I}|\) for \({\mathbb {A}}_{2, \exp }\), and \(1.8 \times R^{-1} |\mathcal {I}|\) for \({\mathbb {A}}_{3, {\mathrm {opt}}}\), where \(R=16 \pi \) is the beam-to-voxel cross-section ratio [see (5.23)].
An isovalue of 12 Gy allows to visualize the locations outside the PTV where the delivered dose is greater than 75% of the target prescription dose, which should be avoided as much as possible even though there is no constraint outside the PTV and the OAR.
Since the dose plan is zero in the OAR, the mean and upper quartile are better indicators of performance than the minimum and standard deviation.
To give an idea, consider for simplicity a spherical tumor of radius r (mm), 1 \(\text {mm}^3\) isotropic voxels, and a depth resolution of 1 mm. In this case \({\mathscr {D}}\) occupies about \(30 \, (r\rho _{\max })^2 |\mathcal {T}| \) bytes of memory, where \(\rho _{\max }\) is the beam radius. So for a 60 mm diameter tumor, a 4 mm beam radius, and 120 GB allocated to \({\mathscr {D}}\), the irradiation geometry is limited to approximately \(1.5 \times 10^5\) trajectories.
Parallel OARs tolerate high doses in small regions provided the mean dose received is small, whereas high doses are harmful to serial OARs even when confined to small regions. Examples of parallel OARs include the lung, kidney and liver; examples of serial OARs include the brain stem, spinal cord and esophagus.
References
Schardt, D., Elsässer, T., Schulz-Ertner, D.: Heavy-ion tumor therapy: physical and radiobiological benefits. Rev. Mod. Phys. 82(1), 383–425 (2010)
Kraft, G., Weber, U.: Tumor therapy with ion beams, pp. 1179–1205. Springer, Handbook of particle detection and imaging (2012)
Loeffler, J., Durante, M.: Charged particle therapy—optimization, challenges and future directions. Nat. Rev. Clin. Oncol. 10(7), 411–424 (2013)
Schlaff, C., Krauze, A., Belard, A., O’Connell, J., Camphausen, K.: Bringing the heavy: carbon ion therapy in the radiobiological and clinical context. Radiat. Oncol. 9, 88 (2014)
Durante, M., Orecchia, R., Loeffler, J.: Charged-particle therapy in cancer: clinical uses and future perspectives. Nat. Rev. Clin. Oncol. 14, 483–495 (2017)
Shepard, D., Ferris, M., Olivera, G., Rockwell Mackie, T.: Optimizing the delivery of radiation therapy to cancer patients. SIAM Rev. 41(4), 721–744 (1999)
Romeijn, H., Ahuja, R., Dempsey, J., Kumar, A., Li, J.: A novel linear programming approach to fluence map optimization for intensity modulated radiation therapy treatment planning. Phys. Med. Biol. 48(21), 3521–3542 (2003)
Romeijn, H., Ahuja, R., Dempsey, J., Kumar, A.: A new linear programming approach to radiation therapy treatment planning problems. Oper. Res. 54(2), 201–216 (2006)
Chen, W., Craft, D., Madden, T., Zhang, K., Kooy, H., Herman, G.: A fast optimization algorithm for multicriteria intensity modulated proton therapy planning. Med. Phys. 37(9), 4938–4945 (2010)
Cao, W., Lim, G., Li, X., Li, Y., Zhu, X., Zhang, X.: Incorporating deliverable monitor unit constraints into spot intensity optimization in intensity-modulated proton therapy treatment planning. Phys. Med. Biol. 58(15), 5113–5125 (2013)
Lomax, A.: Intensity modulation methods for proton radiotherapy. Phys. Med. Biol. 44(1), 185–205 (1999)
Krämer, M., Jäkel, O., Haberer, T., Kraft, G., Schardt, D., Weber, U.: Treatment planning for heavy-ion radiotherapy: physical beam model and dose optimization. Phys. Med. Biol. 45(11), 3299–3317 (2000)
Oelfke, U., Bortfeld, T.: Inverse planning for photon and proton beams. Med. Dosim. 26(2), 113–124 (2001)
Lahanas, M., Schreibmann, E., Baltas, D.: Multiobjective inverse planning for intensity modulated radiotherapy with constraint-free gradient-based optimization algorithms. Phys. Med. Biol. 48(17), 2843–2871 (2003)
Bourhaleb, F., Marchetto, F., Attili, A., Pittà, G., Cirio, R., Donetti, M., Giordanengo, S., Givehchi, N., Iliescu, S., Krengli, M., La Rosa, A., Massai, D., Pecka, A., Pardo, J., Peroni, C.: A treatment planning code for inverse planning and 3D optimization in hadrontherapy. Comput. Biol. Med. 38(9), 990–999 (2008)
Gemmel, A., Hasch, B., Ellerbrock, M., Weyrather, W., Krämer, M.: Biological dose optimization with multiple ion fields. Phys. Med. Biol. 53(23), 6691–7012 (2008)
Unkelbach, J., Bortfeld, T., Martin, B., Soukup, M.: Reducing the sensitivity of IMPT treatment plans to setup errors and range uncertainties via probabilistic treatment planning. Med. Phys. 36(1), 149–163 (2009)
Li, Y., Zhang, X., Mohan, R.: An efficient dose calculation strategy for intensity modulated proton therapy. Phys. Med. Biol. 56(4), N71–N84 (2011)
Inaniwa, T., Kanematsu, N., Furukawa, T., Noda, K.: Optimization algorithm for overlapping-field plans of scanned ion beam therapy with reduced sensitivity to range and setup uncertainties. Phys. Med. Biol. 56(6), 1653–1669 (2011)
Pflugfelder, D., Wilkens, J., Nill, S., Oelfke, U.: A comparison of three optimization algorithms for intensity modulated radiation therapy. Z. Med. Phys. 18(2), 111–119 (2008)
Liu, W., Zhang, X., Li, Y., Mohan, R.: Robust optimization of intensity modulated proton therapy. Med. Phys. 39(2), 1079–1091 (2012)
Romeijn, H., Ahuja, R., Dempsey, J., Kumar, A.: A column generation approach to radiation therapy treatment planning using aperture modulation. SIAM J. Optim. 15(3), 838–862 (2005)
Aleman, D., Glaser, D., Romeijn, H., Dempsey, J.: Interior point algorithms: guaranteed optimality for fluence map optimization in IMRT. Phys. Med. Biol. 55(18), 5467–5482 (2010)
Preciado-Walters, F., Langer, M., Rardin, R., Thai, V.: Column generation for IMRT cancer therapy optimization with implementable segments. Ann. Oper. Res. 148, 65–79 (2006)
Fredriksson, A.: A characterization of robust radiation therapy treatment planning methods—from expected value to worst case optimization. Med. Phys. 39(8), 5169–5181 (2012)
Hartmann, L., Bogner, L.: Investigation of intensity-modulated radiotherapy optimization with gEUD-based objectives by means of simulated annealing. Med. Phys. 35(5), 2041–2049 (2008)
Cotrutz, C., Xing, L.: Segment-based dose optimization using a genetic algorithm. Phys. Med. Biol. 48(18), 2987–2998 (2003)
Li, Y., Yao, J., Yao, D.: Genetic algorithm based deliverable segments optimization for static intensity-modulated radiotherapy. Phys. Med. Biol. 48(20), 3353–3374 (2003)
Ahmad, S., Bergen, S.: A genetic algorithm approach to the inverse problem of treatment planning for intensity-modulated radiotherapy. Biomed. Signal Process. Control 5(3), 189–195 (2010)
Cao, W., Lim, G., Liao, L., Li, Y., Jiang, S., Li, X., Li, H., Suzuki, K., Zhu, X., Gomez, D., Zhang, X.: Proton energy optimization and reduction for intensity-modulated proton therapy. Phys. Med. Biol. 59(21), 6341–6354 (2014)
Cao, W., Lim, G., Lee, A., Li, Y., Liu, W., Zhu, X., Zhang, X.: Uncertainty incorporated beam angle optimization for IMPT treatment planning. Med. Phys. 39(8), 5248–5256 (2012)
Lim, G., Kardar, L., Cao, W.: A hybrid framework for optimizing beam angles in radiation therapy planning. Ann. Oper. Res. 217, 357–383 (2014)
Gu, W., O’Connor, D., Nguyen, D., Yu, V., Ruan, D., Dong, L., Sheng, K.: Integrated beam orientation and scanning-spot optimization in intensity-modulated proton therapy for brain and unilateral head and neck tumors. Med. Phys. 45(4), 1338–1350 (2018)
Ehrgott, M., Güler, Ç., Hamacher, H., Shao, L.: Mathematical optimization in intensity modulated radiation therapy. Ann. Oper. Res. 175, 309–365 (2010)
Craft, D.: Local beam angle optimization with linear programming and gradient search. Phys. Med. Biol. 52(7), N127–N135 (2007)
Aleman, D., Romeijn, H., Dempsey, J.: A response surface approach to beam orientation optimization in intensity-modulated radiation therapy treatment planning. INFORMS J. Comput. 21(1), 62–76 (2009)
Zhang, H., Shi, L., Meyer, R., Nazareth, D., D’Souza, W.: Solving beam-angle selection and dose optimization simultaneously via high-throughput computing. INFORMS J. Comput. 21(3), 427–444 (2009)
Mišić, V., Aleman, D., Sharpe, M.: Neighborhood search approaches to non-coplanar beam orientation optimization for total marrow irradiation using IMRT. Eur. J. Oper. Res. 205(3), 522–527 (2010)
Lee, C.-H., Aleman, D., Sharpe, M.: A set cover approach to fast beam orientation optimization in intensity modulated radiation therapy for total marrow irradiation. Phys. Med. Biol. 56(17), 5679–5695 (2011)
Rocha, H., Dias, J., Ferreira, B., Lopes, M.: Selection of intensity modulated radiation therapy treatment beam directions using radial basis functions within a pattern search methods framework. J. Global Optim. 57(4), 1065–1089 (2013)
Djajaputra, D., Wu, Q., Wu, Y., Mohan, R.: Algorithm and performance of a clinical IMRT beam-angle optimization system. Phys. Med. Biol. 48(19), 3191–3212 (2003)
Aleman, D., Kumar, A., Ahuja, R., Romeijn, H., Dempsey, J.: Neighborhood search approaches to beam orientation optimization in intensity modulated radiation therapy treatment planning. J. Global Optim. 42(4), 587–607 (2008)
Bertsimas, D., Cacchiani, V., Craft, D., Nohadani, O.: A hybrid approach to beam angle optimization in intensity-modulated radiation therapy. Comput. Oper. Res. 40(9), 2187–2197 (2013)
Dias, J., Rocha, H., Ferreira, B., Lopes, M.: Simulated annealing applied to IMRT beam angle optimization: a computational study. Phys. Med. 31(7), 747–756 (2015)
Hou, Q., Wang, J., Chen, Y., Galvin, J.: Beam orientation optimization for IMRT by a hybrid method of the genetic algorithm and the simulated dynamics. Med. Phys. 30(9), 2360–2367 (2003)
Li, Y., Yao, J., Yao, D.: Automatic beam angle selection in IMRT planning using genetic algorithm. Phys. Med. Biol. 49(10), 1915–1932 (2004)
Lei, J., Li, Y.: An approaching genetic algorithm for automatic beam angle selection in IMRT planning. Comput. Methods Programs Biomed. 93(3), 257–265 (2009)
Li, Y., Yao, D., Yao, J., Chen, W.: A particle swarm optimization algorithm for beam angle selection in intensity-modulated radiotherapy planning. Phys. Med. Biol. 50(15), 3491–3514 (2005)
Bangert, M., Ziegenhein, P., Oelfke, U.: Characterizing the combinatorial beam angle selection problem. Phys. Med. Biol. 57(20), 6707–6723 (2012)
Bangert, M., Ziegenhein, P., Oelfke, U.: Comparison of beam angle selection strategies for intracranial IMRT. Med. Phys. 40(1) (2013)
Ehrgott, M., Johnston, R.: Optimization of beam directions in intensity modulated radiation therapy planning. OR Spectr. 25(2), 251–264 (2003)
Lee, E., Fox, T., Crocker, I.: Integer programming applied to intensity-modulated radiation therapy treatment planning. Ann. Oper. Res. 119, 165–181 (2003)
Yang, R., Dai, J., Yang, Y., Hu, Y.: Beam orientation optimization for intensity-modulated radiation therapy using mixed integer programming. Phys. Med. Biol. 51(15), 3653–3666 (2006)
Schreibmann, E., Lahanas, M., Xing, L., Baltas, D.: Multiobjective evolutionary optimization of the number of beams, their orientations and weights for intensity-modulated radiation therapy. Phys. Med. Biol. 49(5), 747–770 (2004)
Fiege, J., McCurdy, B., Potrebko, P., Champion, H., Cull, A.: PARETO: a novel evolutionary optimization approach to multiobjective IMRT planning. Med. Phys. 38(9), 5217–5229 (2011)
Jia, X., Men, C., Lou, Y., Jiang, S.: Beam orientation optimization for intensity modulated radiation therapy using adaptive \(l_{2,1}\)-minimization. Phys. Med. Biol. 56(19), 6205–6222 (2011)
Robini, M., Zhu, Y., Liu, W., Magnin, I.: A stochastic framework for spot-scanning particle therapy. In: Proc. Int. Conf. IEEE Eng. Med. Biol. Soc., pp. 2578–2581, Orlando, FL (2016)
Robini, M.: Theoretically Grounded Acceleration Techniques for Simulated Annealing, Handbook of Optimization, from Classical to Modern Approach, pp. 311–336. Springer (2012)
Robini, M., Reissman, P.-J.: From simulated annealing to stochastic continuation: a new trend in combinatorial optimization. J. Global Optim. 56(1), 185–215 (2013)
Altobelli, E., et al.: Combinatorial optimization in radiotherapy treatment planning. AIMS Med. Sci. 5(3), 204–223 (2018)
Hajek, B.: Cooling schedules for optimal annealing. Math. Oper. Res. 13(2), 311–329 (1988)
Chiang, T.-S., Chow, Y.: On the convergence rate of annealing processes. SIAM J. Control Optim. 26(6), 1455–1470 (1988)
Catoni, O.: Rough large deviation estimates for simulated annealing: application to exponential schedules. Ann. Probab. 20(3), 1109–1146 (1992)
Robini, M., Magnin, I.: Optimization by stochastic continuation. SIAM J. Imaging Sci. 3(4), 1096–1121 (2010)
Weber, U., Kraft, G.: Design and construction of a ripple filter for a smoothed depth dose distribution in conformal particle therapy. Phys. Med. Biol. 44(11), 2765–2775 (1990)
Allison, J., et al.: Geant4 developments and applications. IEEE Trans. Nucl. Sci. 53(1), 270–278 (2006)
Owen, H., Lomax, A., Jolly, S.: Current and future accelerator technologies for charged particle therapy. Nucl. Instrum. Methods Phys. Res. A 809, 96–104 (2016)
Drzymala, R., Mohan, R., Brewster, L., Chu, J., Goitein, M., Harms, W., Urie, M.: Dose-volume histograms. Int. J. Radiat. Oncol. Biol. Phys. 21(1), 71–78 (1991)
Scherrer, A., Küfer, K.-H., Bortfeld, T., Monz, M., Alonso, F.: IMRT planning on adaptive volume structures—a decisive reduction in computational complexity. Phys. Med. Biol. 50(9), 2033–2053 (2005)
Martin, B., Bortfeld, T., Castañon, D.: Accelerating IMRT optimization by voxel sampling. Phys. Med. Biol. 52(24), 7211–7228 (2007)
Catoni, O., Trouvé, A.: Parallel annealing by multiple trials: a mathematical study, Simulated annealing: parallelization techniques, pp. 129–143. Wiley (1992)
Witte, E., Chamberlain, R., Franklin, M.: Parallel simulated annealing using speculative computation. IEEE Trans. Parallel Distrib. Syst. 2(4), 483–494 (1991)
Sohn, A.: Generalized speculative computation of parallel simulated annealing. Ann. Oper. Res. 63, 29–55 (1996)
Karger, C., Peschke, P.: RBE and related modeling in carbon-ion therapy. Phys. Med. Biol. 63(1) (2018)
Thieke, C., Bortfeld, T., Küfer, K.-H.: Characterization of dose distributions through the max and mean dose concept. Acta Oncol. 41(2), 158–161 (2002)
Acknowledgements
This work was supported by the Chinese Basic Research Program under Grant 61671049 and the Chinese Key R&D Program under Grant 2017YFB1400100. We would like to thank Dr. François Smekens for providing us with the cumulative depth-dose profiles of carbon ions in water.
Author information
Authors and Affiliations
Corresponding authors
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
About this article
Cite this article
Robini, M.C., Yang, F. & Zhu, Y. A stochastic approach to full inverse treatment planning for charged-particle therapy. J Glob Optim 77, 853–893 (2020). https://doi.org/10.1007/s10898-020-00902-2
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10898-020-00902-2