pyMKM: An Open-Source Python Package for Microdosimetric Kinetic Model Calculation in Research and Clinical Applications

1. Introduction

Cancer treatment has witnessed remarkable advancements over the years, with radiotherapy standing out as one of the most effective modalities. Among these, hadrontherapy has emerged as a cutting-edge approach, leveraging the unique physical and biological properties of ions to deliver high-dose with exceptional precision. Unlike conventional photon-based radiotherapy, hadrontherapy mitigates damage to healthy tissues while enhancing tumor targeting due to its increased LET, resulting in superior therapeutic outcomes.

The growing interest in hadrontherapy is driven by advancements in accelerator and beam delivery technologies, the development of efficient and comprehensive TPSs, and its increasing clinical adoption. These innovations have consistently demonstrated the efficacy of hadrontherapy, making it a promising alternative to conventional methods. However, realizing its full potential requires the precise assessment of beam quality and the reliable inclusion of the biological effectiveness of different radiation types—key factors that ensure optimal treatment outcomes.

Despite the availability of several models to predict biological effectiveness [1], the MKM [2,3] and its stochastic extension, the SMK model [4,5], have demonstrated superior accuracy in predicting cell survival fractions in vitro, in vivo, and in patients [6]. More recently, the OSMK model [5,7] has been proposed to further extend this framework to hypoxic environments, a critical factor in the treatment of radioresistant tumors. However, their practical application can be challenging due to the complexity and high computational demands of microdosimetric calculations. To address these challenges, we introduce pyMKM, an open-source Python package specifically designed to streamline the generation of microdosimetric tables required by these models. By providing user-friendly scriptable functions for these calculations, pyMKM facilitates the use of MKM, SMK, and OSMK, making them accessible to researchers and clinicians alike.

The development of pyMKM is motivated by the increasing need for centralized and reliable microdosimetric research and clinical platforms, as well as MC codes, such as Geant4 [8], PHITS [9] and FLUKA [10], which have incorporated MKM-related functionalities [4,11,12,13]. The open-source nature of pyMKM promotes collaboration and iterative improvement, enhancing its utility for studies focusing on RBE of ions.

In this paper, we validate pyMKM against published data, including benchmarks from the NIRS, where the models were originally developed as part of their computational framework, and other institutions, to demonstrate its reliability and accuracy. We also outline the key concepts and formulas underpinning MKM, SMK, and OSMK, emphasizing their practical applications within the package.

2. Materials and Methods

This section provides an overview of the MKM and its stochastic extension (SMK model), as well as the OSMK model, focusing on key concepts and formulas necessary for understanding the implementation of the classes within the pyMKM package. These concepts are directly connected to practical examples within pyMKM, including the calculation of survival curves, modeling of specific energy deposition, and the application of both stochastic and oxygen-related corrections for high-LET scenarios.

The structure of this section is twofold. In Sections 2.1–2.4, we consolidate and unify the core formalism that has emerged over the past two decades, gathering in a single place all variables and equations needed for accurate computation. This harmonized presentation facilitates reproducibility and model validation. Detailed discussions of the mathematical formalism and physical-biological principles underpinning these models are provided in the referenced publications and review papers, such as [14]. In Section 2.5, we describe how these theoretical components are mapped onto a modular and test-driven software architecture.

2.1. MKM

The MKM predicts cell survival fractions by modeling the effects of radiation at a microscopic scale, where specific energy deposition in sub-nuclear domains determines the biological outcome. In the context of MKM, domains refer to sub-nuclear microscopic volumes where energy deposition occurs. The nucleus represents the larger cellular structure composed of multiple domains, and its specific energy aggregates the contributions from individual domains to predict the overall survival fraction of the cell. A domain is to be considered dead when a lethal lesion is formed at the site and a cell is to be considered dead when a domain in its nucleus has died [4].

The model assumes that the average number of lethal lesions in the domain follows a linear-quadratic dependence on the specific-energy ($z_d$). The specific energy of the domain is the sum of the single-event specific energy ($z_{d,1}$). Assuming a Poisson distribution of lethal lesions among the irradiated cells, the survival fraction of a cell ($S$) is equal to the probability that the number of lethal lesions is zero in the cell nucleus:

$$S=\exp \left(-L_n\right),$$

(1)

where $L_n$ is the average number of lethal lesions in a nucleus. When the Poisson distribution is assumed for the event number of the domain, $L_n$ is calculated as follows:

$$L_n=-\ln S\left(D\right)={\alpha }_{\mathrm{M}\mathrm{K}\mathrm{M}}D+{\beta }_{\mathrm{M}\mathrm{K}\mathrm{M}}{D}^2,$$

(2a)

$$\begin{array}{l}{\alpha }_{\mathrm{M}\mathrm{K}\mathrm{M}}\equiv {\alpha }_0+{\beta }_0{\overline{z}}_d\\ {\beta }_{\mathrm{M}\mathrm{K}\mathrm{M}}\equiv {\beta }_0\end{array},$$

(2b)

where ${\alpha }_0$ represents the initial slope of the surviving fraction curve at the limit of vanishing LET, $D$ is the absorbed (macroscopic) dose and ${\beta }_0$ is a constant independent of the radiation type [2]. In the case of high LET irradiation, ${\alpha }_0$ and ${\beta }_0$ can be approximated with the experimental ${\alpha }_x$ and ${\beta }_x$ of the reference photon radiation1 [14]. The variable ${\overline{z}}_d$ is the dose-mean specific energy of the domain:

$${\overline{z}}_d={\displaystyle \frac{{\int }_0^{\infty }{z}_{d,1}^2{f}_{d,1}\left(z_{d,1}\right)\mathrm{d}{z}_{d,1}}{{\int }_0^{\infty }{z}_{d,1}{f}_{d,1}\left(z_{d,1}\right)\mathrm{d}{z}_{d,1}}},$$

(3)

where $f_{d,1}\left(z_{d,1}\right)$ is the probability density of $z_{d,1}$ deposited by a single energy-deposition event of the domain [15].

To address the saturation effects observed in high-LET radiation, modifications to the MKM introduced a saturation-corrected specific energy per event ($z_{d,1,\mathrm{s}\mathrm{a}\mathrm{t}}$, also named $z_{d,1}^{\text{'}}$) to account for the overkill effect [2,3]:

$$z_{d,1,\mathrm{s}\mathrm{a}\mathrm{t}}={\displaystyle \frac{z_0^2}{z_{d,1}}}\left[1-\exp \left(-{\displaystyle \frac{z_{d,1}^2}{z_0^2}}\right)\right],$$

(4)

where $z_0$ is the saturation parameter indicating the specific energy above which the correction for the saturation due to the overkill effect becomes significant. The numerical value of $z_0$ is determined based on the assumption of a spherical geometry for the domain as follows:

$$z_0={\displaystyle \frac{{\left(R_n/r_d\right)}^2}{\sqrt{{\beta }_0\left[1+{\left(R_n/r_d\right)}^2\right]}}},$$

(5)

where $R_n$ and $r_d$ are the radii of the cell nucleus and the domain, respectively [3].

Then, ${\overline{z}}_d$ in the linear term of dose in Equation 2b was replaced by the saturation-corrected dose-mean specific energy of the domain (${\overline{z}}_d^{*}$, also referred to as $z_{1\mathrm{D}}^{*}$) [3,4]:

$$-\ln S\left(D\right)={\alpha }_{\mathrm{M}\mathrm{K}\mathrm{M}}^{*}D+{\beta }_{\mathrm{M}\mathrm{K}\mathrm{M}}{D}^2,$$

(6a)

$$\begin{array}{l}{\alpha }_{\mathrm{M}\mathrm{K}\mathrm{M}}^{*}\equiv {\alpha }_0+{\beta }_0{\overline{z}}_d^{*}\\ {\beta }_{\mathrm{M}\mathrm{K}\mathrm{M}}\equiv {\beta }_0\end{array},$$

(6b)

where

$${\overline{z}}_d^{*}={\displaystyle \frac{{\int }_0^{\infty }{z}_{d,1,\mathrm{s}\mathrm{a}\mathrm{t}}{z}_{d,1}{f}_{d,1}\left(z_{d,1}\right)\mathrm{d}{z}_{d,1}}{{\int }_0^{\infty }{z}_{d,1}{f}_{d,1}\left(z_{d,1}\right)\mathrm{d}{z}_{d,1}}}.$$

(7)

Equations 6 correspond to the base implementation of MKM in the pyMKM codebase. Since this modified version, including the saturation correction, is the most commonly adopted, it is not unusual to encounter the expression ${\alpha }_{\mathrm{M}\mathrm{K}\mathrm{M}}\equiv {\alpha }_0+{\beta }_0{\overline{z}}_d^{*}$, directly, without explicitly referencing the transition from Equation 2b to Equation 6b.

It is worth mentioning that some authors, such as [5], define the saturation-corrected specific energy differently, using the following quantity:

$$z_{d,1}^{\text{'}}=z_0\sqrt{1-\exp \left(-{\displaystyle \frac{z_{d,1}^2}{z_0^2}}\right)},$$

(8)

so that Equation 7 becomes

$${\overline{z}}_d^{*}={\displaystyle \frac{{\int }_0^{\infty }{z_{d,1}^{\text{'}}}^2{f}_{d,1}\left(z_{d,1}\right)\mathrm{d}{z}_{d,1}}{{\int }_0^{\infty }{z}_{d,1}{f}_{d,1}\left(z_{d,1}\right)\mathrm{d}{z}_{d,1}}}.$$

(9)

From the perspective of the pyMKM package, adopting the definition in Equation 8 instead the one in Equation 4, offers practical advantages for implementation. Specifically, the use of $z^{\text{'}}=\sqrt{z\cdot z_{\mathrm{s}\mathrm{a}\mathrm{t}}}$ allows for a unified approach to calculating $\overline{z}$ and ${\overline{z}}^{*}$, as their mathematical structures become identical apart from the substitution of $z$ with $z\text{'}$. This unification simplifies the codebase by reducing redundancy and avoids the use of the more complex fraction in the definition of $z_{\mathrm{s}\mathrm{a}\mathrm{t}}$ and the related computational burden of handling division by zero or very small values.

2.2. SMK

The SMK model was developed to address discrepancies between MKM predictions and measured survival fractions for high-LET and high-dose irradiations. While the MKM assumes a constant specific energy for the nucleus ($z_n=D$), this approximation becomes invalid at high LET due to increasing variance in $z_n$ caused by fewer energy deposition events. Consequently, the quadratic coefficient tends to decrease at very high LET, unlike the constant value predicted by the MKM [16].

The SMK incorporates stochastic fluctuations in energy deposition at both the domain ($z_d$) and nucleus ($z_n$) levels [4]. Despite its accuracy [17], the SMK is computationally intensive, requiring significant memory and processing time to account for stochastic variations. To address these challenges, reasonable approximations, as outlined in [5,14], have been introduced to reduce computational demands while maintaining predictive accuracy. With this modified version of the SMK model, the natural logarithm of the survival fraction is calculated as

$$-\ln S\left(D\right)={\alpha }_{\mathrm{S}\mathrm{M}\mathrm{K}}D+{\beta }_{\mathrm{S}\mathrm{M}\mathrm{K}}{D}^2-\ln \left(1+{\gamma }_{\mathrm{S}\mathrm{M}\mathrm{K}}D\right),$$

(10a)

$$\begin{array}{l}{\alpha }_{\mathrm{S}\mathrm{M}\mathrm{K}}\equiv {\alpha }_{\mathrm{M}\mathrm{K}\mathrm{M}}^{*}={\alpha }_0+{\beta }_0{\overline{z}}_d^{*}\\ {\beta }_{\mathrm{S}\mathrm{M}\mathrm{K}}\equiv {\displaystyle \frac{{\overline{z}}_d^{*}}{{\overline{z}}_d}}{\beta }_{\mathrm{M}\mathrm{K}\mathrm{M}}={\displaystyle \frac{{\overline{z}}_d^{*}}{{\overline{z}}_d}}{\beta }_0\\ {\gamma }_{\mathrm{S}\mathrm{M}\mathrm{K}}\equiv {\overline{z}}_n\left[{\displaystyle \frac{1}{2}}{\left({\alpha }_{\mathrm{S}\mathrm{M}\mathrm{K}}+2{\beta }_{\mathrm{S}\mathrm{M}\mathrm{K}}D\right)}^2-{\beta }_{\mathrm{S}\mathrm{M}\mathrm{K}}\right]\end{array},$$

(10b)

where ${\overline{z}}_d$ and ${\overline{z}}_d^{*}$ are the dose-mean and the dose-mean saturation-corrected specific energies absorbed by the domain, and ${\overline{z}}_n$ is the dose-mean specific energy absorbed by the cell nucleus:

$${\overline{z}}_n={\displaystyle \frac{{\int }_0^{\infty }{z}_{n,1}^2{f}_{n,1}\left(z_{n,1}\right)\mathrm{d}{z}_{n,1}}{{\int }_0^{\infty }{z}_{n,1}{f}_{n,1}\left(z_{n,1}\right)\mathrm{d}{z}_{n,1}}},$$

(11)

$f_{n,1}\left(z_{n,1}\right)$ being the single-event probability density of $z_{n,1}$.

As addressed by [4], the saturation parameter $z_0$, introduced in the SMK model when computing ${\overline{z}}_d^{*}$ (Equations 8 and 9), represents the decrease in the number of complex DNA damages per dose, whereas in the MKM, it describes the saturation of cell killing caused by the non-Poisson distribution of lethal lesions among cell nuclei. Consequently, the numerical value of $z_0$ in the SMK model does not necessarily coincides with the one derived from Equation 5 in the MKM. Instead, it is a free parameter of the model.

2.2. OSMK

To extend the SMK formalism to hypoxic conditions2, the OSMK model has been formulated [5]. The model incorporates hypoxia-induced radioresistance by introducing a relative radioresistance factor $\mathcal{R}$, which modulates the production of sublethal lesions as a function of the local oxygen partial pressure $\left[p{O}_2\right]$:

$${\mathcal{R}}_{\left[p{O}_2\right]}={\displaystyle \frac{\left[p{O}_2\right]+K}{\left[p{O}_2\right]+K/{\mathcal{R}}_{\max }}},$$

(12)

where ${\mathcal{R}}_{\max }$ is the maximum $\mathcal{R}$ in the limit of $\left[p{O}_2\right]=0\mathrm{m}\mathrm{m}\mathrm{H}\mathrm{g}$ and $K$ is the oxygen partial pressure at which $\mathcal{R}=({\mathcal{R}}_{\max }+1)/2$.

While lethal lesion formation is assumed to be unaffected by oxygen tension, the production rate of sublethal lesions is reduced by a factor ${\mathcal{R}}^{-1}$, reflecting the higher efficiency of chemical repair mechanisms under low-oxygen conditions.

In this framework the linear term of the survival curve, in the limit of vanishing LET, is decomposed into two components, namely:

$${\alpha }_0={\alpha }_L+{\alpha }_S.$$

(13)

The term ${\alpha }_L$ is associated with lethal lesions formed directly by ionizing radiation; the term ${\alpha }_S$ is associated with sublethal lesions that may evolve into lethal outcomes via biological processes such as interaction or spontaneous conversion. Under hypoxic conditions, only the sublethal component is affected by oxygen tension, while the lethal component remains constant. Analogously, the quadratic coefficient ${\beta }_0$, which scales with the square of the sublethal lesion yield, is reduced more strongly. These substitutions are applied in the SMK model, so that, in the OSMK model, the survival fraction of cells exposed to ionizing radiation under an oxygen partial pressure $\left[p{O}_2\right]$ becomes:

$$-\ln S\left(D\right)={\alpha }_{\mathrm{O}\mathrm{S}\mathrm{M}\mathrm{K}}D+{\beta }_{\mathrm{O}\mathrm{S}\mathrm{M}\mathrm{K}}{D}^2-\ln \left(1+{\gamma }_{\mathrm{O}\mathrm{S}\mathrm{M}\mathrm{K}}D\right),$$

(14a)

$$\begin{array}{l}{\alpha }_{\mathrm{O}\mathrm{S}\mathrm{M}\mathrm{K}}\equiv {\alpha }_L+\left({\alpha }_S/{\mathcal{R}}_{\left[p{O}_2\right]}\right)+{\overline{z}}_d^{*}\left({\beta }_0/{\mathcal{R}}_{\left[p{O}_2\right]}^2\right)\\ {\beta }_{\mathrm{O}\mathrm{S}\mathrm{M}\mathrm{K}}\equiv {\beta }_{\mathrm{S}\mathrm{M}\mathrm{K}}/{\mathcal{R}}_{\left[p{O}_2\right]}^2=\left({\overline{z}}_d^{*}/{\overline{z}}_d\right)\left({\beta }_0/{\mathcal{R}}_{\left[p{O}_2\right]}^2\right)\\ {\gamma }_{\mathrm{O}\mathrm{S}\mathrm{M}\mathrm{K}}\equiv {\overline{z}}_n\left[{\displaystyle \frac{1}{2}}{\left({\alpha }_{\mathrm{O}\mathrm{S}\mathrm{M}\mathrm{K}}+2{\beta }_{\mathrm{O}\mathrm{S}\mathrm{M}\mathrm{K}}D\right)}^2-{\beta }_{\mathrm{O}\mathrm{S}\mathrm{M}\mathrm{K}}\right]\end{array},$$

(14b)

where ${\alpha }_{\mathrm{O}\mathrm{S}\mathrm{M}\mathrm{K}}$ and ${\beta }_{\mathrm{O}\mathrm{S}\mathrm{M}\mathrm{K}}$ appear as ${\alpha }_{\mathrm{S}\mathrm{M}\mathrm{K}}^{\text{'}}$ and ${\beta }_{\mathrm{S}\mathrm{M}\mathrm{K}}^{\text{'}}$, with a prime (‘) symbol used a superscript in the reference literature [5,7].

2.2.1. OSMK Version 2021

In the original formulation [5], ${\mathcal{R}}_{\max }$ is modeled as a function of ${\overline{z}}_d$, thereby incorporating LET-dependent modulation of the OER:

$${\mathcal{R}}_{\max }={\displaystyle \frac{{\left(\raisebox{1ex}{${\overline{z}}_d$}\!\left/ \!\raisebox{-1ex}{$z_{\mathcal{R}}$}\right.\right)}^{{\gamma }_{\mathcal{R}}}+{\mathcal{R}}_m}{{\left(\raisebox{1ex}{${\overline{z}}_d$}\!\left/ \!\raisebox{-1ex}{$z_{\mathcal{R}}$}\right.\right)}^{{\gamma }_{\mathcal{R}}}+1}},$$

(15)

where $z_{\mathcal{R}}$, ${\gamma }_{\mathcal{R}}$ and ${\mathcal{R}}_m$ (i.e. the value of ${\mathcal{R}}_{\max }$ for radiations with very low ${\overline{z}}_d$) are free parameters that describe the physicochemical factors limiting the efficiency of chemical repair mechanisms.

2.2.2. Version 2023

To improve theoretical accuracy in cases of broad LET distributions – such as those found in therapeutic spread-out Bragg peaks – a refined “event-by-event” formulation has been proposed in 2023 [7]. In this version, the hypoxia correction is applied to each individual energy-deposition event rather than to an averaged value. While the same principle of oxygen-dependent reduction in sublethal lesion formation is maintained, key modifications include treating ${\mathcal{R}}_{\max }$ as a cell-specific constant and incorporating direct modulation of the domain radius $r_d$ and saturation parameter $z_0$ with decreasing oxygen pressure. This refinement reflects the hypothesis that hypoxia reduced the effective target volume and diminishes saturation effects even at high LET, leading to:

$$\begin{array}{l}{r}_d^{\text{'}}=r_d/f_{r_d}\\ z_0^{\text{'}}=z_0\cdot f_{z_0}\end{array},$$

(16)

where $f_{r_d}$ and $f_{z_0}$ are two oxygen-dependent modulation factors, functions of the relative radioresistance ${\mathcal{R}}_{\left[p{O}_2\right]}$:

$$\begin{array}{l}{f}_{r_d}\left(\left[p{O}_2\right]\right)=1+\left({\displaystyle \frac{f_{r_d}^{\max }-1}{{\mathcal{R}}_{\max }-1}}\right)\left({\mathcal{R}}_{\left[p{O}_2\right]}-1\right)\\ f_{z_0}\left(\left[p{O}_2\right]\right)={\left[1+\left({\displaystyle \frac{f_{z_0}^{\max }-1}{{\mathcal{R}}_{\max }-1}}\right)\left({\mathcal{R}}_{\left[p{O}_2\right]}-1\right)\right]}^2\end{array},$$

(17)

with $f_{r_d}^{\max }$ and ${f_{z_0}^{\max }}^2$ representing the maximal hypoxia-induced corrections to domain size and saturation energy, respectively.

2.3. Track Structure Models

The calculation of the specific energy per event in the nucleus ($z_{n,1}$) and the domain ($z_{d,1}$​) for MKM applications can be achieved through amorphous track structure models [15] for the single-event specific energy density $d_1={ϵ}_1/\rho$ ($\rho$ being the mass density of the traversed target region with mass $m$ and volume $V$), which provide a numerical approach to theoretically derive the microdosimetric quantities. In this framework, the track structure from an ionizing particle is characterized by radial distributions$d_1\left(\mathit{r}\right)$, which are integrated over the target volume (either the nucleus or domain). The local dose contributions are considered nonzero only within a finite range, extending up to the penumbra radius $R_p$:

$$R_p=\gamma E^{\delta }$$

(18)

depending on the ion’s initial kinetic energy $E$, via two empirical parameters $\gamma$ (which varies based on the track structure model) and a fixed exponent $\delta =1.7$. Consequently, the single-event specific energy $z_1$ is defined as the ratio of the imparted energy ${ϵ}_1$to the mass $m$:

$$z_1={\displaystyle \frac{{ϵ}_1}{m}}={\displaystyle \frac{1}{V}}{\int }_V{d}_1\left(\mathit{r}\right)\mathrm{d}\mathit{r}$$

(19)

Following the computational approach outlined by [3], the sensitive volumes are modeled as cylinders of water ($\rho ={10}^{-12}\mathrm{g}\mathsf{\mu }{\mathrm{m}}^{-1}$) with radius $R$ and height $2R$; the incident ions traverse the sensitive volume with uniform probability; the ion trajectory is assumed to travel parallel to the cylinder’s axis, i.e. $d_1\left(\mathit{r}\right)=d_1\left(r\right)$, with an impact parameter $b$, while any changes in its path or velocity during passage are considered negligible. The single-event specific energy is determined by weighting the deposited local dose according to the overlap between the ion track and the target region, so that Equation 19 becomes:

$$z_1\left(b\right)={\displaystyle \frac{1}{\pi R^2}}{\int }_{R_{\min }}^{R_{\max }}{d}_1\left(r\right)\mathrm{d}\mathcal{A}\left(r,R,b\right)$$

(20)

where $\mathcal{A}\left(r,R,b\right)$ is the intersection area function quantifying this overlap and the integral is evaluated within limits that account for the possible cases of full, partial or no intersection between the ion track and the target region:

$$R_{\min }=\max \left(0,b-R\right),R_{\max }=\min \left(b+R,R_p\right)$$

(21)

The intersection area function $\mathcal{A}$ defines the fraction of the local dose contributing to the region, based on the ion’s relative position, incorporating only the portion of the local dose profile that intersects with the target region:

$$\mathcal{A}\left(r,R,b\right)=\left\{\begin{array}{ll}0,& p\ge r+R\\ \pi \min {\left(r,R\right)}^2,& p\le \left|r-R\right|\\ r^2{\cos }^{-1}\left({\displaystyle \frac{\mathsf{\Delta }+2r^2}{2br}}\right)+R^2{\cos }^{-1}\left({\displaystyle \frac{\mathsf{\Delta }+2R^2}{2bR}}\right)-{\displaystyle \frac{1}{2}}\sqrt{\mathsf{\Sigma }\left(\mathsf{\Sigma }-2r\right)\left(\mathsf{\Sigma }-2R\right)\left(\mathsf{\Sigma }-2b\right)}& \mathrm{o}/\mathrm{w}\end{array}\right.,$$

(22)

where $\mathsf{\Sigma }=r+R+b$ and $\mathsf{\Delta }=b^2-r^2-R^2$.

In pyMKM, two amorphous track structure models were implemented [15]. For both models, it turns out that a fundamental constant can be defined as ${\lambda }_0=1/\left(\pi \rho \right)$:

• Kiefer-Chatterjee model [18], $d_1\left(r\right)=d_{\mathrm{K}\mathrm{C}}\left(r\right)$

$$d_{KC}\left(r\right)=\left\{\begin{array}{ll}{\displaystyle \frac{{\lambda }_0}{R_c^2}}\left[{\mathrm{L}\mathrm{E}\mathrm{T}}_{\mathrm{\infty }}-2\pi \rho K_p\ln \left(R_p/R_c\right)\right],& r\le R_c\\ K_p/r^2,& R_c<r\le R_p\\ 0& r>R_p\end{array}\right.,$$

(23a)

$$R_c=R_0{\beta }_{\mathrm{r}\mathrm{e}\mathrm{l}},K_p=K_0{\left(Z_{\mathrm{e}\mathrm{f}\mathrm{f}}/{\beta }_{\mathrm{r}\mathrm{e}\mathrm{l}}\right)}^2,$$

(23b)

where $R_0=11.6\cdot {10}^{-3}\mathsf{\mu }\mathrm{m}$ is the initial core radius constant parameters, ${\mathrm{L}\mathrm{E}\mathrm{T}}_{\mathrm{\infty }}$ is the unrestricted LET in water, $K_0=1.25\cdot {10}^{-4}\mathrm{G}\mathrm{y}\mathsf{\mu }{\mathrm{m}}^2$, $Z_{\mathrm{e}\mathrm{f}\mathrm{f}}=Z\left[1-\exp \left(-125\cdot {\beta }_{\mathrm{r}\mathrm{e}\mathrm{l}}{Z}^{-{\displaystyle \frac{2}{3}}}\right)\right]$ is the effective charge given by the Barkas expression and ${\beta }_{\mathrm{r}\mathrm{e}\mathrm{l}}={\left(1-{\displaystyle \frac{1}{{\left(1+E/m_0\right)}^2}}\right)}^{1/2}$ is the particle’s velocity relative to the speed of light, with $m_0=931.5\mathrm{M}\mathrm{e}\mathrm{V}/u$ being the atomic rest mass unit. For $d_{\mathrm{K}\mathrm{C}}$, $R_p$ in Equation 18 is calculated with $\gamma =6.16\cdot {10}^{-2}\mathsf{\mu }\mathrm{m}{\left(\mathrm{M}\mathrm{e}\mathrm{V}/u\right)}^{-\delta }$.

• Scholz-Kraft model [19], $d_1\left(r\right)=d_{\mathrm{S}\mathrm{K}}\left(r\right)$

$$d_{\mathrm{S}K}\left(r\right)=\left\{\begin{array}{ll}\lambda \cdot {\mathrm{L}\mathrm{E}\mathrm{T}}_{\mathrm{\infty }}/R_c^2,& r\le R_c\\ \lambda \cdot {\mathrm{L}\mathrm{E}\mathrm{T}}_{\mathrm{\infty }}/r^2,& R_c<r\le R_p\\ 0& r>R_p\end{array}\right.,$$

(24a)

$$\lambda ={\displaystyle \frac{{\lambda }_0}{1+2\ln \left(R_p/R_c\right)}}.$$

(24b)

In pyMKM, different historical parametrization of this model are supported: $R_c={1.0\cdot 10}^{-2}\mathsf{\mu }\mathrm{m}$ (Scholz-Kraft, [19]), $R_c={3.0\cdot 10}^{-4}\mathsf{\mu }\mathrm{m}$ (Elsässer-Sholz, [20]) and $R_c={6.5\cdot 10}^{-2}\mathsf{\mu }\mathrm{m}$ (Friedrich, [21]), with the corresponding coefficient $\gamma$ in Equation 18 set to $\gamma =5.0\cdot {10}^{-2}\mathsf{\mu }\mathrm{m}{\left(\mathrm{M}\mathrm{e}\mathrm{V}/u\right)}^{-\delta }$ (Scholz-Kraft, [19]) or $\gamma =6.2\cdot {10}^{-2}\mathsf{\mu }\mathrm{m}{\left(\mathrm{M}\mathrm{e}\mathrm{V}/u\right)}^{-\delta }$ (Elsässer-Sholz, [20] or Friedrich, [21]).

Independently of the chosen track structure, pyMKM applies a lower bound on the penumbra radius by enforcing $R_p\ge R_c$. This condition prevents the penumbra from becoming smaller than the high-dose core and ensures numerical stability of the radial dose integration. In practice, the cutoff becomes effective only at very low kinetic energies, typically below ~0.1 MeV/u, while in the clinically relevant range the calculated values of $R_p$ remain unaffected.

2.4. Calculation of Dose-Averged Quantities

Dose-mean saturation-corrected and uncorrected quantities are computed by weighting the corresponding specific energy per event according to its probability of occurrence, $f_1$, as described in Equations 3, 9 and 11, which share the same functional form.

When considering a monoenergetic ion beam, the specific energy per event follows the deterministic relationship $z_1=z_1\left(b\right)$, dictated by the track structure model (Equation 19). In this case, the dose-averaged integrals can be rewritten in terms of the impact parameter $b$:

$$\overline{z}={\displaystyle \frac{{\int }_0^{\infty }{z}_1^2{f}_1\left(b\right)\mathrm{d}b}{{\int }_0^{\infty }{z}_1{f}_1\left(b\right)\mathrm{d}b}}$$

(25)

The justification for this transformation lies in the uniform transverse distribution of fluence assumed. In a parallel ion beam, ion trajectories are uniformly distributed in the transverse plane, meaning that the probability of an ion having an impact parameter in the range $\left[b,b+\mathrm{d}b\right]$ is proportional to the area of the annular region at radius $b$, i.e. $f_1\left(b\right)\mathrm{d}b\propto 2\pi b\mathrm{d}b$. As a result, the integral for $\overline{z}$ can be rewritten in terms of $b$ [3], incorporating the weighting factor $b$ that naturally arises from the cylindrical geometry of the ion fluence distribution:

$$\overline{z}={\displaystyle \frac{{\int }_0^{b_{\max }}{z}_1^2\left(b\right)b\mathrm{d}b}{{\int }_0^{b_{\max }}{z}_1\left(b\right)b\mathrm{d}b}},$$

(26)

where

$$b_{\max }=R_p+R.$$

(27)

Here, $b_{\max }$ represents the maximum impact for which a non-negligible energy deposition occurs in the target region. It is given by the sum of the penumbra radius of the track $\left(R_p\right)$ and the radius of the sensitive volume $\left(R\right)$.

2.5. Software Implementation

The pyMKM package is implemented in Python using a modular architecture (Figure 1), systematically mirroring the theoretical structure of MKM, SMK, and OSMK models presented in Sections 2.1–2.4. Its design separates concerns between radiobiological modeling, particle transport and geometry, microdosimetric computation, and I/O data. The following describes the internal structure and the main classes available in each module.

2.5.1. Microdosimetric Table Generation (pymkm.mktable)

The generation of microdosimetric tables is managed through the MKTable class, located in pymkm.mktable.core. It encapsulates the model version, an MKTableParameters object, and a reference to the corresponding StoppingPowerTableSet. The MKTableParameters class stores all user-defined inputs, such as the radiobiological model, domain and nucleus radii, number of sampling points, and flags for stochastic and oxygen effect corrections. The compute submodule integrates track structure calculations with biologically relevant quantities such as dose-averaged specific energies and saturation-corrected values. It communicates with the particle_track and specific_energy modules from pymkm.physics, and supports parallelization through the parallel utility.

2.5.2. Survival Table Computation (pymkm.sftable)

Survival curve generation is handled by the SFTable class within pymkm.sftable.core, which uses inputs defined in a SFTableParameters instance. It requires a pre-computed MKTable object and computes the natural logarithm of the survival fraction over a configurable dose grid. The table structure supports multi-ion mixtures and stochastic corrections (SMK) by reusing modules from mktable.compute. The class hierarchy ensures that SFTableParameters maintains a reference to the corresponding MKTable, facilitating consistent parameter propagation. The results can be visualized using the plot module.

2.5.3. Track Structure and Energy Deposition (pymkm.physics)

Two fundamental components enable microdosimetric calculations:

• ParticleTrack: Encapsulates the radial dose distribution model, LET, particle energy, and stopping power information. It provides access to initial local dose and penumbra radius calculations, serving as the primary physical abstraction for energy deposition events;
• SpecificEnergy: Computes single-event and dose-averaged specific energy distributions, with or without saturation correction. It interfaces with ParticleTrack to perform event-based computations, supporting MKM, SMK, and OSMK formalisms.

These components also interact with the geometry_tools module for evaluating dose contributions as a function of the impact parameter and radial geometry (Equation 22).

2.5.4. Oxygen Effect Modeling (pymkm.biology.oxygen_effect)

To enable modeling under hypoxic conditions, the oxygen_effect module provides implementations of the 2021 and 2023 versions of the OSMK model. It supports LET-dependent modulation of $\alpha$ and$\beta$, as well as event-level correction of domain radius and saturation parameters in the latest formulation. These are automatically applied when relevant flags are activated in MKTableParameters or SFTableParameters.

2.5.5. Stopping Power and Data Management (pymkm.io and pymkm.data)

The StoppingPowerTable and StoppingPowerTableSet classes manage the import, interpolation, and organization of stopping power datasets. These are extracted from MSTAR (https://nds.iaea.org/stopping-legacy/MstarWWW/MSTARInstr.html), FLUKA or Geant4 simulations and are accessible from multiple submodules across the package.

StoppingPowerTable supports I/O in *.json, *.txt, and internal formats, and provides energy and stopping power values for arbitrary ion species. StoppingPowerTableSet aggregates multiple tables and performs batch operations such as interpolation, filtering, and energy range extraction. It also acts as the key data source within MKTable. The data.defaults submodules provide preloaded datasets for reproducibility and rapid prototyping.

2.5.6. Utilities and Shared Components (pymkm.utils)

• GeometryTools offers static methods for sampling radial positions, computing overlapping areas, and generating resolution-optimized radii.
• Interpolation is used across modules to resample LET and energy distributions, with support for linear and log-log interpolation.
• Parallel streamlines multiprocessing routines, ensuring efficient computation of event-specific microdosimetric quantities.

2.6. Testing, Validation, and Continuous Integration

The pyMKM package is developed following principles of test-driven development and reproducible science. A comprehensive suite of automated unit and integration tests covers all core functionalities. These tests can be executed locally using standard pytest invocation.

In addition to the testing infrastructure, the package includes a collection of example scripts and model-specific validation scripts used to benchmark pyMKM against published data. These scripts serve both as usage demonstrations and as tools for verifying consistency with theoretical and experimental results from the literature. In total, the validation suite comprises over 150 reference datasets spanning different ion species, energies, LET values, and experimental conditions, with comparisons against published benchmarks from NIRS.

The package has been validated only within the energy range of 0.1 MeV/u to 1000 MeV/u, which covers the scope of currently available reference data and encompasses the energies of clinical relevance in ion beam therapy.

Each validation script reproduces key reference figures or tables from the literature and verifies that pyMKM reproduces the expected trends in local dose, specific energy, MKM/SMK/OSMK survival predictions, and table generation workflows. All comparisons are accompanied by quantitative error metrics—such as absolute or logarithmic deviations—adapted to the physical scale and domain of the reference data.

To ensure long-term reliability, all tests and benchmarks are also continuously executed in clean environments via GitHub Actions. This continuous integration setup automatically verifies the integrity of the codebase upon each commit or pull request, ensuring platform-independent reproducibility and robustness of the computational results. Tests are executed in clean environments using multiple Python versions, guaranteeing cross-platform stability and maintainability of the package.

3. Results

The reference values for assessing the accuracy of pyMKM were retrieved using a semi-automatic data digitization pipeline based on the GSYS graphical data grabber (https://www.jcprg.org/gsys/2.4/gsys-e.html). As a result, all comparisons inherently include uncertainties associated with graphical data extraction, which may contribute to small systematic or localized deviations when evaluating the numerical agreement.

While selected cases are presented in this section to illustrate representative outcomes, the validation relies on quantitative metrics computed over a benchmark suite comprising approximately 150 datasets. Complete results—including all validation figures and corresponding numerical error tables—are available in the official GitHub repository at https://github.com/BeppeMagro/pymkm

3.1. Validation of Local Dose and Specific Energy Distributions

The validation campaign covered ion beams from hydrogen to carbon ($Z=1\mathrm{t}\mathrm{o}6$), with LET values ranging from ~10 to over 500 keV/μm. For each test case, pyMKM outputs were compared against reference values via interpolation and error analysis.

For initial local dose (Figure 2), quantitative agreement with the digitized reference data was consistently good across all tested ion species and LET values. The mean logarithmic error was 0.20 ± 0.06 log₁₀ units, with RMS deviations of 0.40 ± 0.09 and maximum localized deviations up to 1.48 log₁₀ units. Despite these outliers, the overall relative error remained low, with an average SMAPE of 2.8% (range: 1.7–3.4%).

For specific energy distributions (Figure 3), accuracy improved further. The mean logarithmic error was reduced to 0.083 ± 0.034, with RMS errors of 0.17 ± 0.08 and maxima typically below 0.6. The corresponding SMAPE averaged 8.4%, although with a broader spread (0.7–17%).

Finally, for saturation-corrected specific energy, the results were comparable to the uncorrected case. The mean log error remained at 0.083 ± 0.023, with RMS values of 0.17 ± 0.05 and maxima around 0.55. The SMAPE averaged 7.9% (range: 1.8–17.2%), confirming the numerical stability of the overkill formalism even in high-LET regimes. The similarity of the error statistics between the corrected and uncorrected distributions indicates that the saturation modification is propagated consistently through the integration routines, without introducing additional systematic bias.

To illustrate the overall agreement, Figure 2, Figure 3 and Figure 4 show three representative cases: (a) radial dose distribution comparison for a high-LET carbon ion using both track structure models; (b) specific energy as a function of impact parameter for different sensitive region sizes; and (c) saturation-corrected profiles across three domain/nucleus configurations. These examples were selected to reflect both the physical range of configurations tested and the structural accuracy of the implementation.

In all cases, pyMKM predictions (dashed lines) show excellent agreement with the digitized reference data (solid lines), demonstrating the numerical accuracy of the radial dose integration and microdosimetric kernel generation modules.

3.2. Validation of Microdosimetric Table Generation

To validate the MKTable implementation in pyMKM, we benchmarked the computed dose-averaged domain-specific energy $\left({\overline{z}}_d\right)$, the saturation-corrected dose-averaged domain-specific energy $\left({\overline{z}}_d^{*}\right)$, and the dose-averaged nucleus-specific energy $\left({\overline{z}}_n\right)$ against digitized reference values extracted from published literature.

For the saturation-corrected dose-averaged domain-specific energy, the agreement with reference values was excellent across the entire energy range. The mean absolute error was 0.29 ± 0.05 Gy, with maximum deviations of 1.6 ± 0.7 Gy. Relative differences remained consistently below 2%, with an average SMAPE of 1.4% (range: 1.2–1.7%).

The stochastic implementation was evaluated across multiple combinations of ion species and geometrical configurations. The overall mean absolute error was 4.0 ± 7.5 Gy, with maximum deviations occasionally exceeding 20 Gy, while the average SMAPE remained low at 2.5% ± 1.2% (range: 1.1–5.4%). A more detailed inspection of the three fundamental microdosimetric quantities highlighted the following:

• ${\overline{z}}_d^{*}$: showed the best agreement, with mean absolute errors in the range 0.2–0.5 Gy, maximum deviations below 3 Gy, and SMAPE values of 1–2% across all species;

• ${\overline{z}}_d$: absolute deviations increased with Z, with mean errors rising from ~1 Gy (Z = 1) to ~24 Gy (Z = 10) and maxima up to ~88 Gy. Despite these large absolute numbers, relative differences remained low, with SMAPE values between 1.8% and 4.6%;

• ${\overline{z}}_n$: absolute values were much smaller (<1 Gy), so errors also remained on the order of 0.001–0.03 Gy. In this narrow scale, relative errors appeared slightly larger (SMAPE 2–5%), but overall agreement with reference data remained excellent.

When stratified by atomic number, light ions (H, He) showed somewhat higher fluctuations, while heavier ions (C, Ne) maintained smoother agreement.

Finally, the validation of the survival-relevant $D₁₀$ dose confirmed the robustness of pyMKM. The mean absolute error was 0.014 ± 0.005 Gy, with maxima below 0.08 Gy. Relative differences were consistently <0.6%.

Figure 5 and Figure 6 exemplify these comparisons and demonstrate that both the magnitude and shape of microdosimetric trends are faithfully reproduced by pyMKM. This supports the internal consistency of track structure models, saturation correction logic, and table computation workflows.

3.3. Validation of Survival Fraction Tables

The accuracy of the SFTable module in pyMKM was evaluated across a comprehensive set of survival curves generated under various conditions, including different ion types, LET values, and radiobiological models (MKM, SMK, OSMK). A total of 96 cases were benchmarked against digitized reference data, covering normoxic and hypoxic conditions.

For MKM, validation included 18 survival curves across HSG and V79 cells (Figure 7). In HSG cells, agreement with reference data was excellent, with mean log errors typically below 0.03 log₁₀ units, maxima ≤0.10, and SMAPE values consistently <1.2%. In contrast, V79 cells exhibited systematically larger deviations, with mean log errors up to 0.12, maxima around 0.27, and SMAPE values between 2–3.5%. Overall, errors remained below 0.1 log₁₀ units and ~2% SMAPE for carbon and neon in HSG cells, while deviations in V79 cells were higher but still within acceptable limits for survival modeling.

The stochastic extension improved accuracy across all conditions (Figure 8). Mean log errors were reduced to ~0.016 ± 0.008 log₁₀ units, with maximum deviations never exceeding 0.06. SMAPE values remained below 2% in all cases and often <1%. The improvement was most pronounced at high LET. Importantly, the largest discrepancies observed in V79 cells under MKM were significantly reduced in SMK, with maximum deviations dropping below 0.06 log units and relative differences ≤2%. HSG and CHO results were uniformly accurate.

Validation under hypoxic conditions was carried out using two complementary approaches:

• Direct survival curves at $\left[p{O}_2\right]=0\mathrm{m}\mathrm{m}\mathrm{H}\mathrm{g}$ (Figure 9). Across 42 benchmarked cases, the mean log error was 0.024 ± 0.013 and SMAPE 1.9% ± 0.8%, with maxima <0.16 log₁₀ units. The agreement held across ions and cell types, although H and He again showed slightly higher fluctuations in V79 cells.
• OER vs LET at 10% survival (Figure 10). For 18 benchmark cases, the mean absolute error was 0.025 ± 0.018, maxima <0.07, and SMAPE values below 1%. The characteristic reduction of OER with increasing LET was faithfully reproduced for all cell lines.

Across 96 benchmarked cases, pyMKM reproduced survival data with excellent accuracy: MKM within ~2%, SMK within ~1%, and OSMK within ~2% for hypoxic conditions, with OER predictions accurate to better than 1%.

Figure 7, Figure 8, Figure 9 and Figure 10 provide representative examples of the validation campaign. In particular, Figure 7 shows MKM survival curves for carbon ions across increasing LET; Figure 8 presents SMK-based curves for neon ions; Figure 9 highlights survival under hypoxia using the OSMK model; and Figure 10 illustrates the agreement in OER behavior. Although the figures illustrate representative cases, the main validation rests on the consistent quantitative agreement observed across all tested scenarios.

3.4. Software Validation and Test Coverage

In addition to model-specific validation against reference biological data, the internal consistency and structural robustness of pyMKM were assessed through a complete suite of unit and integration tests.

At the time of this release, the testing infrastructure included 273 individual test cases. All tests were executed using the pytest 8.3.4 framework on a Windows 10 platform with Python 3.9.23 and completed successfully with zero failures. Test coverage analysis reported 100% statement coverage across ~1500 lines of code, with no uncovered branches or functions. This level of coverage reflects not only the completeness of the testing infrastructure, but also the modular design of pyMKM, which facilitates precise targeting of each computational component.

To ensure long-term reproducibility and platform independence, all tests were also executed via continuous integration on GitHub Actions, as part of an automated testing workflow. The CI pipeline verified that all test cases pass consistently across multiple Python versions, specifically 3.9, 3.10, 3.11, and 3.12, with execution times below one minute per environment. This cross-version compatibility further strengthens confidence in the package’s maintainability and future extensibility.

These results confirm that pyMKM is not only numerically accurate in reproducing radiobiological outcomes, but also structurally sound and robust. The automated test suite provides continuous verification of code correctness and guards against unintended side effects during development or model extension.

4. Discussion

The results presented in this work validate the full computational pipeline of pyMKM across a broad spectrum of radiobiological conditions, using survival data digitized from multiple published sources. The benchmark set included different ion species, LET values, cell lines, and radiobiological models (MKM, SMK, OSMK), covering both normoxic and hypoxic environments. Microdosimetric tables and survival curves were compared against reference values using multiple metrics, confirming the high numerical accuracy of the implementation in reproducing model predictions.

All stopping power data used in this validation were derived from the MSTAR dataset, selected for practical reasons due to its reduced sampling density and faster computation. However, the pyMKM framework also supports stopping power inputs from FLUKA and Geant4 MC codes, which are available and validated internally. The consistency of results across the full radiobiological modeling chain—from stopping power to survival fraction—demonstrated the reliability of the track structure integration routines and the robustness of the numerical pipeline even when relying on a simplified input.

The validation of pyMKM against a broad set of datasets confirms the correct numerical implementation of MKM, SMK, and OSMK, but also highlights systematic features that clarify the origin of certain deviations. For local dose and specific energy distributions, the apparently inconsistent behavior of error metrics is explained by their different sensitivities. Radial dose profiles showed mean logarithmic errors around 0.20 log₁₀, higher than for specific energy (~0.08), yet their SMAPE values were smaller (~3% versus ~8%). This inversion arises because the logarithmic metric accentuates discrepancies in the low-value tails of the radial profile, which span several orders of magnitude, while SMAPE normalizes differences against the local magnitude. Similarly, specific energy and saturation-corrected distributions occasionally exhibited SMAPE values up to 17%, even when curves overlapped visually, because small denominators in regions of very low absolute values inflate relative deviations. Importantly, these regions often correspond to asymptotic approaches to zero at low LET or in the far radial tails, where the semi-automatic digitization of reference curves is particularly prone to error. In such cases, even modest uncertainties in data extraction may generate disproportionately large deviations, reinforcing the need to interpret quantitative metrics alongside visual inspection.

For microdosimetric table generation, the behavior of the error metrics strongly reflects the magnitude of the underlying quantities. Domain-specific energies showed the largest absolute deviations, increasing from ~1 Gy for protons to more than 20 Gy for neon and argon, with maxima exceeding 80 Gy at Z = 10. Nevertheless, relative errors (SMAPE 1–4%) remained consistently low, confirming that the agreement is excellent when normalized to the large scale of these values. Saturation-corrected domain quantities were even more stable, with absolute deviations typically below 0.5 Gy and SMAPE values close to 1–2%, indicating that the overkill correction is propagated consistently through the calculation. At the nucleus scale, the opposite situation arises: absolute errors were extremely small (0.001–0.03 Gy), but SMAPE values fluctuated between 2% and 5%. This apparent deterioration is explained by the small magnitude of the reference values themselves, which remain below ~1.3 Gy even for the heaviest ion tested. In this low-value regime, tiny numerical differences are disproportionately amplified by relative metrics, although they are negligible from a biological standpoint. Thus, while the metrics appear less stable for nucleus-specific quantities, the underlying agreement remains excellent and consistent with expectations.

Importantly, pyMKM also reproduced the characteristic U-shaped dependence of the dose at 10% survival as a function of LET. This result is particularly relevant because D₁₀ represents a macroscopic and clinically meaningful endpoint, widely used in radiobiology to compare ion species and to benchmark models against experimental survival data. The fact that pyMKM correctly reproduces this quantity confirms not only the internal consistency of the saturation correction but also the reliability of the package in capturing trends that directly impact clinical interpretation.

Survival fraction validation further confirmed the reliability of pyMKM but also exposed a systematic anomaly limited to V79 cells. The results reported here already account for a uniform scaling factor of 1.3 applied to the reference LET, without which agreement with survival curves could not be achieved. The remarkable point is that this correction seems to have propagated across all ions in the benchmark V79 datasets once it was first introduced for helium. A plausible explanation is that early measurements with He-3 beams were later reported on an He-4 LET scale, leading to a ~30% shift in LET values. When these data were subsequently used as reference, the mis-scaled LET axis appears to have been applied consistently to the entire set of V79 survival curves, not just helium. This would explain why the anomaly is exclusive to V79 but affects all ions. Importantly, the systematic nature of the shift points to a dataset-level inconsistency rather than a limitation of pyMKM’s implementation. Beyond this peculiarity, the model comparisons behaved as expected: the SMK extension performed accuracy at high LET, while OSMK reproduced both hypoxic survival curves and OER–LET relationships, confirming correct incorporation of oxygen-dependent lesion modulation.

Finally, the assessment of software validation and test coverage demonstrated that pyMKM is not only numerically accurate in reproducing radiobiological outcomes, but also structurally sound and robust. The automated test suite, combined with continuous integration across Python versions, ensures continuous verification of code correctness and guards against unintended side effects during development or model extension. This level of reproducibility and platform independence provides an essential foundation for extending the package to future model formulations and clinical applications.

5. Conclusions

We have introduced pyMKM, a fully open-source Python package for the computation of microdosimetric quantities and radiobiological survival predictions based on MKM, SMK, and OSMK models. The package was validated against multiple reference datasets across different ion species, LET values, and oxygenation levels, demonstrating high numerical accuracy and model consistency.

All core modules achieved 100% test coverage and passed continuous integration checks across multiple Python versions, ensuring long-term stability and reproducibility. These results confirm pyMKM as a reliable and extensible framework for microdosimetric modeling in both research and clinical contexts.