0031-9155/70/24/245027

Objective. This study aims to comprehensively compare the PENHAN and FLUKA Monte Carlo codes for low-energy alpha particle transport and small-scale dosimetry using alpha-emitting radionuclides, and to assess their suitability for such applications. Approach. Two studies were performed through Monte Carlo simulations. First, monoenergetic alpha particles ( MeV) were distributed in a micrometric water sphere and the dose deposition within it was calculated. Second, a simplified spherical cell model with uniformly distributed alpha-emitting radionuclides was used to compute S-values. PENHAN and FLUKA results were compared, and potential sources of discrepancy between them were analyzed. In addition, both codes were benchmarked against MIRDcell, an analytical tool widely used for dosimetric calculations in Targeted radionuclide therapy. Main results. In the monoenergetic study, the primary source of discrepancy between PENHAN and FLUKA was the stopping powers used for alpha particles. When the same stopping powers were employed, both codes yielded statistically compatible results, except at 3.0 and 3.5 MeV, where FLUKA showed an anomalous behavior. In the cell model, variations were below 3% but not negligible even when using identical stopping powers, suggesting an additional source of discrepancy: differences in the radionuclide emission spectra, particularly in the electron component. In both studies, PENHAN and FLUKA results were overall in good agreement with those from MIRDcell. Significance. This study demonstrates, for the first time, the suitability of PENHAN for low-energy alpha transport and small-scale dosimetry with alpha emitters, provided that accurate stopping powers are employed. It also supports the reliability of FLUKA in these scenarios and shows that both codes yield compatible results when using consistent stopping power datasets and radionuclide emission spectra. This work further highlights the importance of validating Monte Carlo codes in medical physics to ensure the reliability and reproducibility of their results.

Export citation and abstractBibTeXRIS

Targeted radionuclide therapy (TRT) is a form of radiotherapy in which a radioactive isotope is distributed into the patient’s body through the administration of a radiopharmaceutical, which accumulates in the tumor area. As a result, the radioisotope emissions deposit their energy locally within the tumor, producing lethal damage to cancerous cells while minimizing irradiation to surrounding healthy tissue.

This technique traditionally uses β− electron emitters, such as , and (Salih et al 2022). However, radionuclides that emit low-energy alpha particles (typically in the range of 3–10 MeV) are gaining attention due to their advantageous radiobiological properties. Alpha particles at these energies have shorter ranges in tissue: approximately in contrast with of β− electrons, whose maximum energies typically reach . Additionally, they have a significantly higher linear energy transfer (LET), with values in the range of , compared to the characteristic of β− radiation.

These features imply that alpha radiation produces a more localized dose deposition at the tumor site, sparing the surrounding healthy tissue, and also inflicts more effective damage to cancer cells. Consequently, alpha emitters have emerged as promising candidates for TRT, and several works have already been published, both reviewing the current state of Targeted Alpha Therapy (Eychenne et al 2021, Jang et al 2023, Miederer et al 2024) and describing specific clinical strategies involving different alpha-emitting radionuclides, such as (Albertsson et al 2023), (Hooijman et al 2024, Rahman et al 2024) and (Kokov et al 2022, Scaffidi-Muta and Abell 2025), the latter being a generator of the alpha emitter .

Besides alpha and β− emitters, Auger electron emitters such as , , and have also been investigated in recent years for TRT applications. Auger electrons exhibit even shorter ranges in matter than alpha particles—from nanometers to a few micrometers—and a higher LET than β− electrons—approximately (Ku et al 2019). In addition, these radionuclides enable combined therapeutic and diagnostic use, as Auger electrons are responsible for cellular damage, while the additional photons they emit allow monitoring through external detectors. Several studies have explored the potential of these emitters, often referred to as theranostic agents (Falzone et al 2015, McMillan et al 2015, Kusumoto et al 2020).

When investigating TRT techniques, it is essential to perform small-scale dosimetry to assess dose and energy deposition at the cellular level (Tronchin et al 2022). The reason is that the radiopharmaceuticals used accumulate within or near the target cells, and, especially for alpha and Auger electron emitters, the dose is deposited close to the decay site due to their short ranges. Understanding dose distribution at these scales is therefore critical for accurately evaluating treatment efficacy and potential side effects.

To carry out cellular dosimetry calculations, both theoretical (analytical) models and computational methods are employed, since experimental techniques are limited at the small scales involved. Analytical approaches offer faster computation times but rely on approximations that reduce accuracy; this makes them suitable for clinical applications, where rapid estimations are required, but less appropriate for research, where time constraints are less critical and higher precision is preferred. In contrast, computational methods—specifically Monte Carlo simulations—provide a more detailed description of all radiation-matter interaction processes, accounting for effects such as energy straggling, which are often neglected in analytical models, and allow for the study of complex geometries. Consequently, these techniques have become a powerful and widely used tool in small-scale dosimetry studies (Chatzipapas et al 2020, Levrague et al 2025).

Taking into account the important role that Monte Carlo simulations play in medical physics research, it is essential to validate both the simulation algorithms and the underlying physics models used to describe radiation-matter interactions. This is particularly important in the case of alpha particle transport, for which significantly less research has been conducted compared to electrons, photons, and protons, which are more commonly used in clinical treatments (Mairani et al 2022). In cellular dosimetry studies, where the spatial scales are so small that experimental methods are not applicable, the most effective approach to validate Monte Carlo methods is to compare the results obtained using different codes and ensure their consistency.

Several Monte Carlo codes are available for the simulation of radiation transport, with some being more suitable than others depending on the specific application. In the field of medical physics, the most commonly used general-purpose Monte Carlo codes are (Andreo et al 2017): GEANT4 (Agostinelli et al 2020), FLUKA (Bohlen et al 2014, Battistoni et al 2015, 2016), PENHAN (Salvat and Quesada 2020, 2021, Salvat and Heredia 2024, Salvat et al 2025)—an extension of PENELOPE (Salvat 2023) that includes the transport of neutrons, protons, and alpha particles -, and PHITS (Sato et al 2024). All of these codes allow for the transport of electrons, positrons, photons, and heavy charged particles over the energy ranges relevant to medical physics applications. In this work we use two of them: PENHAN and FLUKA. The previous version of PENHAN, PENELOPE, has been extensively used for dosimetric studies involving emitters used in TRT (Stewart et al 2002, Bernal and Liendo 2009, Hocine et al 2014, Falzone et al 2015, Mendes et al 2021, Valente et al 2025). However, since alpha particle transport was only recently incorporated into PENHAN (in the January 2024 release), there is currently a lack of studies in the literature making use of this new feature, which highlights the relevance of the present work. In contrast, FLUKA’s transport of hadrons, including alpha particles, is well established, as demonstrated by studies available in the literature (Sommerer et al 2006, Tessonnier et al 2017). Nevertheless, its application to small-scale dosimetry has not been extensively explored, with only a few published works addressing this topic (Botta et al 2017, Scarinci et al 2025). As reference data to compare the results obtained with PENHAN and FLUKA, we use the MIRDcell software (version v.14.6) (Vaziri et al 2014, Katugampola et al 2022, 2024), which employs an analytical approach to compute dosimetric quantities (Goddu et al 1994). This formalism is recommended by the Medical Internal Radiation Dose (MIRD) Committee and is considered a standard dosimetric methodology for TRT.

The aim of this work is to evaluate the suitability of the PENHAN and FLUKA codes for performing small-scale dosimetry with alpha-emitting radionuclides, as well as to carry out a thorough comparison between them in order to identify the sources of discrepancy leading to differences in their dosimetric calculations. Results obtained with both codes are benchmarked against those provided by the MIRDcell software.

2.1. Monte Carlo codes

In this section, we briefly summarize the main features of the two Monte Carlo codes employed in this work.

2.1.1. PENHAN

PENHAN (Salvat and Quesada 2020, 2021, Salvat and Heredia 2024, Salvat et al 2025) is a Monte Carlo code, developed by Francesc Salvat at the University of Barcelona (Spain), for the simulation of the coupled transport of electrons, positrons, photons, protons, neutrons, and alpha particles in diverse materials for energies ranging from a few hundred eV to 1 GeV. It is an extension of the previously released simulation code PENELOPE (Salvat 2023), which only included electron, positron, and photon transport. In this work, we used the 2024 version of PENHAN.

To calculate the particle trajectories, the algorithm utilizes an interaction model, mainly consisting of a set of differential cross-sections for the involved interaction mechanisms, obtained from both numerical databases and analytical models. The main interaction processes that are taken into account are: photoelectric effect, Rayleigh scattering, Compton scattering and electron–positron pair production for photons; elastic scattering, inelastic collisions, bremsstrahlung emission and positron annihilation for electrons and positrons; elastic collisions, inelastic collisions and nuclear reactions induced by proton impact for heavy particles (protons, neutrons and alphas).

PENHAN is a class-II simulation code, meaning that it uses a mixed approach, where hard interactions are simulated in a detailed mode, while soft interactions are treated using a condensed history method, based on approximate multiple scattering theory. The algorithm classifies interactions as hard or soft events, depending on whether the energy loss and/or angular deflection of the particle being tracked exceeds (hard) or remains below (soft) certain cutoff values.

These cutoff values governing the mixed simulations are controlled by the input parameters , , C1 and C2 regarding electrons, positrons, and photons, and by , and regarding hadrons. Additionally, the input parameters define the transport energy cutoffs for different particle types, i.e. below these thresholds, particle transport is stopped and the remaining energy is deposited locally.

A particularly notable feature of PENHAN is its ability to perform fully detailed simulations by adjusting the input parameters. In this mode, all interactions are treated as hard events; therefore, all particle tracks are explicitly followed, avoiding the approximations and averaging inherent to condensed history algorithms. As a result, the simulation outcomes are more precise, provided that the underlying physical models are properly validated. However, detailed simulations have the drawback of significantly increased computational time and should therefore be reserved for cases where such a high level of detail is necessary.

2.1.2. FLUKA

FLUKA (Bohlen et al 2014, Battistoni et al 2015) is a Monte Carlo code that allows the simulation of the transport of approximately 60 different particle types, including photons, electrons, and positrons from 100 eV–1 keV to thousands of TeV, hadrons with energies up to 20 TeV, neutrons down to thermal energies and heavy ions, such as alpha particles, with energies starting from a few keV. It has been developed by CERN in Geneva (Switzerland) in collaboration with INFN in Italy. In this work, we used FLUKA version 4–4.1.

FLUKA employs a class-I algorithm, meaning it performs condensed history simulations, where particle trajectories are described in steps that encompass multiple interactions. The effects of these interactions are averaged along each step using an algorithm based on Moliére’s Multiple Elastic Scattering theory.

To conduct our simulations, we used the FLAIR Graphical User Interface (version 3.3–1) (Donadon et al 2024), which facilitates the preparation of input files and execution of simulations without requiring extensive programming knowledge. With this tool, we easily defined the necessary CARDS to configure our simulations by setting transport energy cutoffs or energy cutoffs for electron and bremsstrahlung production, selecting the precision of the physical models used, and configuring the scorers required to obtain the physical quantities under investigation.

2.2. Simulation setup2.2.1. Configuration parameters

In the simulations performed with PENHAN, we adopted a mixed approach, with the simulation parameters summarized in columns 3–4 of table 1. The transport parameters were set according to the recommended values provided in the code manual, while the absorption energies were chosen as the minimum values for which the underlying physical models remained valid.

Table 1. Simulation parameters used in PENHAN (columns 3–4) and FLUKA (columns 5–7) for tracking electrons, positrons, photons, and charged hadrons. Italicized values in parentheses in the PENHAN columns indicate the modifications applied to perform fully detailed simulations.

Particle typeParameter descriptionPENHANFLUKAParameterValueCARDParameterValuePhotonTransport energy cutoffEABS(2)EMFCUT with SDUM = blankWHAT(2)Energy threshold for photon production——EMFCUT with SDUM = PROD-CUTWHAT(2)Electron/ positronTransport energy cutoffEABS(1) and EABS(3)EMFCUT with SDUM = blankWHAT(1)Energy threshold for electron and positron production——EMFCUT with SDUM = PROD-CUTWHAT(1)Mixed simulation parametersC10.05 (0)———C20.05 (0)———WCC———WCH (eV)———Charged hadronProton transport energy cutoffEABSH(1)PART-THResWHAT(1)WHAT(2)PROTONAlpha particle transport energy cutoffEABSH(2)PART-THResWHAT(1)WHAT(2)4-HELIUMEnergy threshold for delta ray production——DELTARAYWHAT(1)Mixed simulation parametersC1H0.001 (0)———C2H0.05 (0)———WCCH———AllSelection of default parameters for precision simulations——DEFAULTS with SDUM = PRECISIO——

Additionally, given the small distances and low energies involved in the studied scenarios, we also conducted simulations using a detailed approach to evaluate whether it could provide relevant insights into the study. To perform these detailed simulations, we set the parameters C1, C2, and to 0, and equal to the electron transport energy cutoff, and equal to the photon transport energy cutoff, as shown by the italicized values in parentheses in table 1. All other parameters were kept the same as in the mixed approach.

For FLUKA, we used the PRECISIOn model, which applies the recommended physics settings to our simulations. Additionally, we set the absorption energies to the minimum values allowed by the code, and adjusted the remaining parameters to match—as closely as possible—those used in PENHAN mixed simulations. The input CARDS used, along with their corresponding parameters, are detailed in columns 5–7 of table 1.

In Monte Carlo simulations, the number of primary particles simulated is crucial since it affects the precision of the results. Simulating more particles improves the accuracy but also increases the computational time required. In our simulations, we used between and 106 primary particles, which ensured uncertainties below for the study involving monoenergetic alpha sources and below for the one involving alpha-emitting radionuclides. All uncertainties reported in this work correspond to a coverage factor of k = 3. No variance reduction techniques were applied in the simulations performed.

All simulations were run on a Core i7-1270P 12th Gen (12 cores, 16 threads), except for those performed with PENHAN using the detailed approach, which were executed on an AMD Epyc ROME 7662 (64 cores, 128 threads) due to their higher computational cost.

2.2.2. Stopping power values

In Monte Carlo simulations, the underlying physics models play a crucial role, as their accuracy determines the reliability of the results.

For hadron transport in matter, several physical interactions can occur, including Coulomb scattering, ionization (i.e. inelastic collisions with atomic electrons), and nuclear interactions such as elastic scattering with atomic nuclei or nuclear reactions that produce other species. Among these, ionization and nuclear reactions are the main processes affecting dose deposition.

In TRT scenarios, the low energies involved make nuclear reactions negligible, and therefore, the main mechanism impacting our dose calculations is ionization. In Monte Carlo simulations, this process is implemented using the electronic stopping power, which quantifies the energy loss of the particle as it traverses the medium. Consequently, the stopping power values used for alpha particles in water in the simulations—and how they are calculated—have a direct impact on the dose calculations.

In the case of PENHAN, stopping powers are calculated using the extended Sternheimer–Liljequist generalized oscillator strength model. However, the code also allows users to input their own stopping power values, enabling the use of externally derived data. This feature ensures that more accurate values can be incorporated when available. Accordingly, in our study, in addition to the PENHAN default values, we also employed a set of stopping powers generated by the SBETHE software (Salvat and Andreo 2023). SBETHE, distributed with the PENHAN code, applies the corrected Bethe formula for energies above specific thresholds (0.75 MeV for protons and 5 MeV for alphas) and empirical relations derived from IAEA experimental data for lower energies. The output is provided on the same energy grid as PENHAN, making the replacement process straightforward.

For FLUKA simulations, we used the default stopping powers provided by the code, which are calculated based on the Bethe–Bloch formalism with several corrections (Battistoni et al 2016), and are consistent, by construction, with those reported in NIST databases PSTAR and ASTAR (Berger et al 2005).

Figure 1 presents the stopping power values for alpha particles in water provided by FLUKA, SBETHE, and PENHAN (default values) over the energy range relevant to TRT scenarios. It also includes experimental data from various studies (Palmer and Akhavan-Rezayat 1978, Akhavan-Rezayat and Palmer 1980, Thwaites 1981, Haque et al 1985), extracted from the IAEA stopping power database (Montanari et al 2024). Both SBETHE and FLUKA are in good agreement with the experimental values in the 1–10 MeV energy range. In contrast, the PENHAN default stopping powers, which clearly deviate from both SBETHE and FLUKA, show significant discrepancies with experimental data at energies below 2 MeV. In this work, we analyzed how these discrepancies affect dose deposition.

Figure 1. Mass stopping powers for low-energy alpha particles in water provided by FLUKA (purple), SBETHE (green), and PENHAN default (orange). Lines are included as a visual guide. Experimental values obtained from the IAEA stopping power database (Montanari et al 2024), based on the studies by Palmer and Akhavan-Rezayat (1978), Akhavan-Rezayat and Palmer (1980), Haque et al (1985), and Thwaites (1981), are also shown (various shades of gray). The inset is included to show the complete low-energy behavior for the PENHAN curve, which is truncated in the main figure, in comparison with FLUKA and SBETHE.

Download figure:

Standard image High-resolution image 2.2.3. Simulation scenarios

We considered two simplified geometrical scenarios. In the first case, corresponding to the study with monoenergetic alpha particle sources, a water sphere of radius was placed inside a larger water sphere of radius, serving as an approximation of an infinite water medium. Alpha particles with energies ranging from 3–—typical in TRT treatments—were uniformly distributed in the inner sphere.

In the second scenario, which involved alpha-emitting radionuclides relevant to TRT, a simplified model of a spherical cell was used. This model consisted of two concentric water spheres with radii of and , representing the nucleus and the cytoplasm, respectively. As in the previous case, the cell model was embedded within a radius water sphere. Two source distributions were considered: one with the radionuclide uniformly distributed within the nucleus, and another within the cytoplasm. These two configurations allowed for the calculation of the self-dose to the cell nucleus when the radiopharmaceutical is intracellularly taken up—a common situation in TRT. The compartment in which uptake occurs—nucleus or cytoplasm—depends on the specific chemical properties of the radiopharmaceutical. Other source configurations are also possible, such as a uniform distribution throughout the entire cell (nucleus plus cytoplasm), or localization on the cell membrane surface (Sgouros et al 2021), although these were not considered in the present work.

Liquid water was selected as the medium in all cases, as its properties closely resemble those of soft biological tissues. In the simulations, the pre-defined water material available in each code was employed; both codes used a mean excitation energy of .

2.3. S-value calculation

In small-scale dosimetry studies related to TRT, a key objective is to assess dose deposition at the cellular level when radionuclide uptake occurs within the target cell.

A common approach to this analysis is that proposed by the MIRD committee, involving the S-value. This quantity, denoted as , represents the mean absorbed dose in a target region T per radioactive decay in a source region S. According to the MIRD formalism, the S-value is defined as (Bolch et al 2009):

where ni is the number of particles of type i emitted per decay, Ei is the mean energy of radiation of type i, is the fraction of energy emitted from S as radiation of type i that is absorbed in T and is the mass of the target region.

This quantity is particularly relevant, as it allows estimation of the mean absorbed dose within the target region over a given period of time, provided the cumulated activity inside the source region —i.e. the total number of decays occurring in the source during that period—is known:

In our study, S-values were derived from simulation results through simple calculations. The simulations provided energy deposition per primary particle simulated, and the absorbed dose—and consequently the S-value—was obtained by dividing the deposited energy by the target mass.

The calculations were performed using the simplified cellular model described in the previous section. S-values were determined at the cell nucleus when the radionuclide was uniformly distributed within the nucleus, , and when it was distributed within the cytoplasm, .

2.4. Simulation of radionuclides

We analyzed dose deposition by three radionuclides considered promising candidates for targeted alpha therapy: