Kinetic-diffusion-rotation algorithm for dose estimation in electron beam therapy

Klaas Willems; Vince Maes; Zhirui Tang; Giovanni Samaey

doi:10.3934/krm.2025018

Article Contents

2026, Volume 19: 160-181. Doi: 10.3934/krm.2025018

This volume Previous Article From kinetic mixtures to compressible two-phase flow: A BGK-type model and rigorous derivation Next Article A novel linear transport model with distinct scattering mechanisms for direction and speed

Kinetic-diffusion-rotation algorithm for dose estimation in electron beam therapy

1.
KU Leuven, Belgium
2.
Rheinland-Pfälzische Technische Universität Kaiserslautern-Landau, Germany

^*Corresponding author: Klaas Willems
^*Corresponding author: Klaas Willems

Received: March 2025

Revised: August 2025

Early access: September 25, 2025

Published online: September 25, 2025

Abstract / Introduction Full Text(HTML) Figure(8) / Table(3) Related Papers Cited by

Abstract

Monte Carlo methods are state-of-the-art when it comes to dosimetric computations in radiotherapy. However, the execution time of these methods suffers in high-collisional regimes. We address this problem by introducing a kinetic-diffusion particle tracing scheme. This algorithm, first proposed in the context of neutral transport in fusion energy, relies on the explicit simulation of the exact kinetic motion in low-collisional regimes and dynamically switches to an approximate random walk in high-collisional regimes. The random walk corresponds to an advection-diffusion process that preserves the first two moments (mean and variance) of the kinetic motion. We derive an analytic formula for the mean kinetic motion and discuss the addition of a multiple scattering distribution to the algorithm. In contrast to neutral transport, the electron beam therapy setting does not readily admit to an analytical expression for the variance of the kinetic motion, and we therefore resort to the use of a lookup table. We test the algorithm for dosimetric computations in electron beam therapy on a 2D CT scan of a lung patient. Using a simple particle model, our Python implementation of the algorithm is nearly 33 times faster than an equivalent kinetic simulation at the cost of a small modeling error.

Keywords:

Mathematics Subject Classification: Primary: 35Q20, 60G50and65C05; Secondary: 92-10.

Citation:

Full Text(HTML)

Figure 1. Angular deflections at scattering events [29]

Download: Full-size image PowerPoint slide

Figure 2. Mean kinetic motion (left) and the dependence of the mean polar scattering angle on the energy (right)

Download: Full-size image PowerPoint slide

Figure 3. Comparison of experimentally obtained mean kinetic motion (10.000 particles) with the theoretical result (13) as a function of the energy. Particles are set to travel a distance $ \Delta s = 0.1 $

Download: Full-size image PowerPoint slide

Figure 4. Lookup table for the variance in the x, y and z direction generated with a kinetic particle tracing algorithm. Specific simulation parameters are summarized in table 1. The error bars represent one standard deviation computed using the batch means method with 10 batches [7]

Download: Full-size image PowerPoint slide

Figure 5. Polar multiple scattering distribution $ \theta_{MS} $ for a particle with energy $ E_{max} = 3.117 $ MeV, that has travelled $ \Delta s = 0.1 $ centimetres in a medium with density $ \rho = 1.0 \; g/cm^3 $. The distribution is fitted with a log-normal distribution

Download: Full-size image PowerPoint slide

Figure 6. Speed-up of KDR compared to purely kinetic simulation for a particle with scattering characteristics of a 2.61 MeV electron

Download: Full-size image PowerPoint slide

Figure 7. 2D CT scan of a lung and the computed dose distribution using three different particle tracing algorithms. The borders of the colored regions are isodose lines for $ 10^{-8}, 10^0, 10^2, 10^3 $, and $ 10^5 \frac{\text{MeV}}{\text{cm}^2} $. All relevant simulation parameters are summarised in table 2

Download: Full-size image PowerPoint slide

Figure 8. Pointwise relative error of the dose distributions

Download: Full-size image PowerPoint slide

Table 1. Parameters used for the generation of the lookup table for the variance of the kinetic motion

$ E_{min} $	$ E_{max} $	$ \rho_{min} $	$ \rho_{max} $	$ \Delta s_{min} $	$ \Delta s_{max} $
0.5 MeV	21 MeV	$ 0.05 \; g/cm^3 $	$ 1.85 \; g/cm^3 $	$ 1\times10^{-4} \; cm $	$ 1 \; cm $

	Amount of bins for $E$, $\rho$ and $\Delta s$			Amount of particles
	16			10⁴

| Show Table

DownLoad: CSV

Table 2. Simulation parameters for the lung test case with KDR. The initial energy of the particle is given by a normal distribution with mean $ E_{max} $ and standard deviation $ \sigma_E $. The domain has the shape of a square with size $ d $ in the y-z direction and is infinite in the x direction. The initial position of the particle is given by (0, $ d $, $ z_i $), where $ z_i $ is sampled from a Gaussian distribution with mean $ d/2 $ and standard deviation $ \sigma_z $. For each simulation algorithm, $ N $ particles are simulated. The KDR stepsize is $ \Delta s $. The initial velocity of the particles is constrained to the y-z plane where the angle with the negative y-axis is given by a von Mises distribution with center zero and dispersion $ \kappa_\theta $

$ E_{max} $	$ \sigma_E $	$ d $	$ \sigma_z $	$ \kappa_\theta $	$ \Delta s $	$ N $
21 MeV	$ \frac{1}{100} $	14.5 $ cm $	$ \frac{1}{50} $	10000	0.0725	100.000

| Show Table

DownLoad: CSV

Table 3. Timings results for the lung test case obtained on the Flemisch supercomputer using 16 MPI processes. Results are given in the hh:mm:ss format

Analog	KDR	KDR MS
06:12:51	00:09:33	00:11:20

| Show Table

DownLoad: CSV

Related Papers

Cited by

References

[1]	A. Al-Beteri and D. Raeside, Designing, benchmarking and applying a Monte Carlo electron transport code, Computer Methods and Programs in Biomedicine, 39 (1983), 147-167. doi: 10.1016/0169-2607(93)90019-H.
[2]	J. K. Blitzstein and J. Hwang, Introduction to Probability, 2^nd edition, CRC Press, Boca Raton, 2019.
[3]	I. Boyd and T.Deschenes, Hybrid Particle-Continuum Numerical Methods for Aerospace Applications, May 2012.
[4]	C. Cercignani, The Boltzmann Equation and Its Applications, No. 67 in Applied mathematical sciences. Springer, New York, NY Heidelberg, 1988. doi: 10.1007/978-1-4612-1039-9_2.
[5]	G. X. Ding, D. M. Duggan, C. W. Coffey, P. Shokrani and J. E. Cygler, First macro Monte Carlo based commercial dose calculation module for electron beam treatment planning—new issues for clinical consideration, Physics in Medicine and Biology, 51 (2006), 2781-2799. doi: 10.1088/0031-9155/51/11/007.
[6]	F. C. P. du Plessis, C. A. Willemse, M. G. Lotter and L. Goedhals, The indirect use of CT numbers to establish material properties needed for Monte Carlo calculation of dose distributions in patients, Medical Physics, 25 (1998), 1195-1201. doi: 10.1118/1.598297.
[7]	J. M. Flegal and G. L. Jones, Batch means and spectral variance estimators in Markov chain Monte Carlo, The Annals of Statistics, 38 (2010), 1034-1070. doi: 10.1214/09-AOS735.
[8]	C. K. Garrett and C. D. Hauck, A comparison of moment closures for linear kinetic transport equations: The line source benchmark, Transport Theory and Statistical Physics, 42 (2013), 203-235. doi: 10.1080/00411450.2014.910226.
[9]	S. Hissoiny, B. Ozell, H. Bouchard and P. Després, GPUMCD: {A} new GPU-oriented Monte Carlo dose calculation platform, Medical Physics, 38 (2011), 754-764. doi: 10.1118/1.3539725.
[10]	K. Jabbari, Review of fast Monte Carlo codes for dose calculation in radiation therapy treatment planning, Journal of Medical Signals and Sensors, 1 (2011), 73.
[11]	A. Jablonski, C. J. Powell and S. Tanuma, Monte Carlo strategies for simulations of electron backscattering from surfaces, Surface and Interface Analysis, 37 (2005), 861-874. doi: 10.1002/sia.2104.
[12]	D. Jacqmin, SU-{D}-218-06: Acceleration of Optical Photon Monte Carlo Simulations Using the Macro Monte Carlo Method, Medical Physics, 39 (2012), 3623-3623. doi: 10.1118/1.4734709.
[13]	I. Kawrakow, M. Fippel and K. Friedrich, 3D electron dose calculation using a Voxel based Monte Carlo algorithm (VMC), Medical Physics, 23 (1996), 445-457. doi: 10.1118/1.597673.
[14]	I. Kawrakow, E. Mainegra-Hing, D. W. O. Rogers, F. Tessier and B. R. B. Walters, The EGSnrc Code System: Monte Carlo Simulation of Electron and Photon Transport.,
[15]	P. J. Keall and P. W. Hoban, Super-Monte Carlo: {A} 3-{D} electron beam dose calculation algorithm, Medical Physics, 23 (1996), 2023-2034. doi: 10.1118/1.597842.
[16]	P. J. Keall and P. W. Hoban, Superposition dose calculation incorporating Monte Carlo generated electron track kernels, Medical Physics, 23 (1996), 479-485. doi: 10.1118/1.597679.
[17]	K. Kisling, L. Zhang, H. Simonds, N. Fakie, J. Yang, R. McCarroll, P. Balter, H. Burger, O. Bogler, R. Howell, K. Schmeler, M. Mejia, B. M. Beadle, A. Jhingran and L. Court, Fully automatic treatment planning for external-beam radiation therapy of locally advanced cervical cancer: {A} tool for low-resource clinics, Journal of Global Oncology, 5 (2019), 1-9. doi: 10.1200/JGO.18.00107.
[18]	J. Kusch and P. Stammer, A robust collision source method for rank adaptive dynamical low-rank approximation in radiation therapy, ESAIM Math. Model. Numer. Anal., 57 (2023), 865-891.
[19]	K. Küpper, Models, Numerical Methods and Uncertainty Quantification for Radiation Therapy., PhD thesis, Aachen University, Nov. 2016.
[20]	E. W. Larsen and J. B. Keller, Asymptotic solution of neutron transport problems for small mean free paths, Journal of Mathematical Physics, 15 (1974), 75-81. doi: 10.1063/1.1666510.
[21]	E. Loevbak, Multilevel and Adjoint Monte Carlo Methods for Plasma Edge Neutral Particle Models, PhD thesis, KU Leuven, Jan. 2023.
[22]	E. Løvbak and G. Samaey, Accelerated simulation of Boltzmann-BGK equations near the diffusive limit with asymptotic-preserving multilevel Monte Carlo, SIAM Journal on Scientific Computing, 45 (2023), A1862-A1889. doi: 10.1137/22M1498498.
[23]	V. Maes, W. Dekeyser, J. Koellermeier, M. Baelmans and G. Samaey, Hilbert expansion based fluid models for kinetic equations describing neutral particles in the plasma edge of a fusion device, Physics of Plasmas, 30 (2023), 063907. doi: 10.1063/5.0146158.
[24]	B. Mortier, Advanced Monte Carlo Simulation and Estimation for Kinetic Neutral Particles in the Plasma Edge of Fusion Reactors, PhD thesis, KU Leuven, 2020.
[25]	B. Mortier, M. Baelmans and G. Samaey, A kinetic-diffusion asymptotic-preserving monte carlo algorithm for the boltzmann-BGK model in the diffusive scaling, SIAM Journal on Scientific Computing, 44 (2022), A720-A744. doi: 10.1137/20M1381526.
[26]	B. Mortier, V. Maes and G. Samaey, Estimation as a post‐processing step for random walk approximations of the Boltzmann‐BGK model, Contributions to Plasma Physics, 62 (2022). doi: 10.1002/ctpp.202100197.
[27]	H. Neuenschwander and E. J. Born, A macro Monte Carlo method for electron beam dose calculations, Physics in Medicine and Biology, 37 (1992), 107-125. doi: 10.1088/0031-9155/37/1/007.
[28]	H. Neuenschwander, T. R. Mackie and P. J. Reckwerdt, MMC-a high-performance Monte Carlo code for electron beam treatment planning, Physics in Medicine and Biology, 40 (1995), 543-574. doi: 10.1088/0031-9155/40/4/005.
[29]	Nuclear Energy Agency, PENELOPE 2018: {A Code System for Monte Carlo Simulation of Electron and Photon transport: Workshop Proceedings, Barcelona, Spain, 28 January – 1 February 2019}, PENELOPE: {A} code system for Monte Carlo simulation of electron and photon transport. OECD, Sept. 2019.
[30]	E. Olbrant, Models and Numerical Methods for Time- and Energy-Dependent Particle Transport, PhD thesis, RWTH Aachen University, 2012.
[31]	H. G. Othmer and T. Hillen, The diffusion limit of transport equations derived from velocity-jump processes, SIAM Journal on Applied Mathematics, 61 (2000), 751-775. doi: 10.1137/S0036139999358167.
[32]	N. Reynaert, S. Van Der Marck, D. Schaart, W. Van Der Zee, C. Van Vliet-Vroegindeweij, M. Tomsej, J. Jansen, B. Heijmen, M. Coghe and C. De Wagter, Monte Carlo treatment planning for photon and electron beams, Radiation Physics and Chemistry, 76 (2007), 643-686. doi: 10.1016/j.radphyschem.2006.05.015.
[33]	D. J. Rhee, A. Jhingran, K. Kisling, C. Cardenas, H. Simonds and L. Court, Automated radiation treatment planning for cervical cancer, Seminars in Radiation Oncology, 30 (2020), 340-347. doi: 10.1016/j.semradonc.2020.05.006.
[34]	S. M. Seltzer, Electron-photon Monte Carlo calculations: The ETRAN code, International Journal of Radiation Applications and Instrumentation. Part A. Applied Radiation and Isotopes, 42 (1991), 917-941. doi: 10.1016/0883-2889(91)90050-B.
[35]	O. N. Vassiliev, Monte Carlo Methods for Radiation Transport, Biological and Medical Physics, Biomedical Engineering. Springer International Publishing, Cham, 2017.
[36]	N. W. R. and N. Yoshihito, The EGS4 Code System: Solution of gamma-ray and electron transport problems.,
[37]	M. B. Wilk and R. Gnanadesikan, Probability Plotting Methods for the Analysis of Data, Biometrika, 55 (1968), 1. doi: 10.2307/2334448.
[38]	K. Willems, ElectronTransportCode: Repository containing implementations of kinetic, kinetic-diffusion and kinetic-diffusion-rotation particle tracing algorithms for radiation therapy, 2023, https://github.com/KlaasWillems/electronTransportCode.
[39]	K. Willems, Particle, Fluid and Hybrid Numerical Methods for Radiation Therapy, Master's thesis, KU Leuven, Leuven, Belgium, 2023.