Integrating microdosimetric in vitro RBE models for particle therapy into TOPAS MC using the MicrOdosimetry-based modeliNg for RBE ASsessment (MONAS) tool

Objective. In this paper, we present MONAS (MicrOdosimetry-based modelliNg for relative biological effectiveness (RBE) ASsessment) toolkit. MONAS is a TOPAS Monte Carlo extension, that combines simulations of microdosimetric distributions with radiobiological microdosimetry-based models for predicting cell survival curves and dose-dependent RBE. Approach. MONAS expands TOPAS microdosimetric extension, by including novel specific energy scorers to calculate the single- and multi-event specific energy microdosimetric distributions at different micrometer scales. These spectra are used as physical input to three different formulations of the microdosimetric kinetic m odel, and to the generalized stochastic microdosimetric model (GSM2), to predict dose-dependent cell survival fraction and RBE. MONAS predictions are then validated against experimental microdosimetric spectra and in vitro survival fraction data. To show the MONAS features, we present two different applications of the code: (i) the depth-RBE curve calculation from a passively scattered proton SOBP and monoenergetic 12C-ion beam by using experimentally validated spectra as physical input, and (ii) the calculation of the 3D RBE distribution on a real head and neck patient geometry treated with protons. Main results. MONAS can estimate dose-dependent RBE and cell survival curves from experimentally validated microdosimetric spectra with four clinically relevant radiobiological models. From the radiobiological characterization of a proton SOBP and 12C fields, we observe the well-known trend of increasing RBE values at the distal edge of the radiation field. The 3D RBE map calculated confirmed the trend observed in the analysis of the SOBP, with the highest RBE values found in the distal edge of the target. Significance. MONAS extension offers a comprehensive microdosimetry-based framework for assessing the biological effects of particle radiation in both research and clinical environments, pushing closer the experimental physics-based description to the biological damage assessment, contributing to bridging the gap between a microdosimetric description of the radiation field and its application in proton therapy treatment with variable RBE.


Introduction
Proton therapy is now widely recognized as an advanced form of radiation therapy compared to the conventional use of photons for treating a steadily increasing number of types of cancers, especially deep-seated, radioresistant, and hypoxic [Loeffler andDurante, 2013, Tambas et al., 2022].The advantages of ions over photons are mainly attributed to the localized energy deposition at the end-of-range of the particle, known as the Bragg peak, resulting in a highly conformal dose distribution and normal tissue sparing [Durante et al., 2017].In addition to the physical advantages, ion beam therapy is characterized by larger biological effectiveness.The reason for this is the higher ionization density and more severe damage to cellular DNA (e.g.double-strand breaks and clustered damage) than photon radiation [Scholz et al., 2001].The superior biological effectiveness of ions is quantified by the Relative Biological Effectiveness (RBE), that is, the ratio between the X-ray (reference) and ion dose causing a same biological effect [Jäkel et al., 2016].In particle therapy, in general, to obtain a homogeneous biological effect in the tumor volume, the RBE of ions is included in the treatment plan optimization as a multiplication factor to the absorbed physical dose.In proton treatment planning, an RBE value equal to 1.1 is adopted for both tumor and normal tissue, namely, protons are considered 10% more effective than photons.Despite the clinical practice assumes a spatially invariant RBE for protons, pre-clinical in vitro and in vivo studies have demonstrated that a constant RBE is an oversimplification due to its dependence on numerous parameters (dose, dose rate, cell line, biological endpoint, radiation quality in the specific voxel, etc.) [Paganetti andGoitein, 2000, Paganetti, 2014], with the RBE being significantly above 1.1 in the distal region, [Missiaggia et al., 2020].For heavier ions like Carbon and Helium, the variations in RBE along the beam penetration depth are so significant that a fixed RBE value cannot be deemed appropriate, and current treatment plans are calculated accounting for a variable RBE, [Inaniwa andKanematsu, 2018, Mairani et al., 2022].
Therefore, for optimal treatment outcomes that effectively balance the targeting of the tumor and the minimization of the damage to the surrounding healthy tissue, it is crucial to accurately estimate the RBE at any point of the irradiated field.The first step in achieving such a challenging goal is to characterize the radiation field, both in terms of macroscopic absorbed energy and the microscopic local pattern of energy deposition.Microdosimetry [Zaider et al., 1996] has proven over the years to be an extremely powerful tool to accomplish such a task.Microdosimetry is a branch of physics that studies the energy deposition of particles at a scale in the order of a few microns, which is the scale of a cell nucleus, believed to be the most sensitive target to radiation-induced cell killing due to the presence of DNA.At the micron scale, singleparticle energy deposition is characterized by large fluctuations due to the inherently stochastic nature of particle interaction, and, therefore, microdosimetry characterizes the radiation field in terms of probability distributions of energy.
By characterizing the energy depositions at the micron scale, microdosimetry provides an ideal tool to link radiation to its biological effects directly.For this reason, many radiobiological models rely on microdosimetry principles, among which the Microdosimetric Kinetic Model (MKM) is the most prominent and widespread in particle therapy [Hawkins, 1994, Hawkins, 1996, Inaniwa and Kanematsu, 2018].
The MKM is a mechanistic model that predicts the cell survival fraction of irradiated cells based on microdosimetric average values and estimates the resulting RBE [Zaider et al., 1996].In particular, the MKM predicts the logarithm of the cell survival fraction of irradiated cells as a linear-quadratic (LQ) function of the imparted dose, [McMahon, 2018], log S(D) = −αD − βD 2 , (1) with α and β two radiobiological parameters that depend both on the biological tissue and on the specific ionizing radiation.Despite the MKM has displayed notable consistency with experimental data obtained both in vitro and in vivo [Mein et al., 2020], over time, numerous successive adaptations of the MKM have been developed in the literature, [Kase et al., 2006, Kase et al., 2007, Inaniwa et al., 2010, Sato et al., 2006, Bellinzona et al., 2021a], primarily aimed to address limitations of the original model in specific scenarios where its underlying assumptions were unsuitable.Currently, the MKM represents the standard model used to calculate the RBE in several carbon-ion therapy centers, [Mein et al., 2020], and, along with the Local Effect Model (LEM), the MKM is one of the only two models currently used in clinics for this purpose.Further, the MKM has been used as the reference RBE model in the recently first treated patient with helium, [Mairani et al., 2022].The MKM's success and widespread use highlight the importance of microdosimetry in accurately predicting the biological effects of radiation and optimizing treatment planning for patients.
Recently, the Generalized Stochastic Microdosimetric Model (GSM 2 ) has been developed, [Cordoni et al., 2021], which is a theoretically grounded mechanistic radiobiological model able to include several spatiotemporal stochastic effects inherent to the formation and repair of radiation-induced DNA damage, [Cordoni et al., 2022b, Cordoni et al., 2022a, Missiaggia et al., 2022].GSM 2 is a fully probabilistic model that overcomes one of the main assumptions shared by most existing radiobiological models including the MKM, that is the fact that the distribution of the number of damages induced by radiation on DNA is Poissonian.In doing so, GSM 2 describes the time evolution of probability distribution of radiation-induced DNA damage rather than focusing on average values as done in the MKM.
Although mechanistic models based on the microdosimetric description of radiation field quality have been shown as an accurate tool for predicting RBE in ion therapy, experimental microdosimetric spectra are still challenging to measure in the daily clinical practice even using commercial detectors, e.g.tissue equivalent proportional counter (TEPC) or silicon-oninsulator (SOI) detector [Bradley et al., 2001, Kase et al., 2006, Rosenfeld, 2016, Missiaggia et al., 2020, Conte et al., 2020, Missiaggia et al., 2021, Missiaggia et al., 2023].Therefore, numerical algorithms, such as Monte Carlo (MC) particle simulation toolkits, have been demonstrated as a valuable alternative for the microdosimetric characterization of the radiation field [Baratto-Roldán et al., 2021, Zhu et al., 2019, Missiaggia et al., 2020, Missiaggia et al., 2023].Nevertheless, MC simulations for microdosimetric spectra are exceedingly timeconsuming, which has prevented their integration into clinical treatment planning systems.This limitation has hindered the utilization of valuable microdosimetric insights within the clinical practice.Hence, to expedite computational processes various numerical approximations have been incorporated into MKM models, [Inaniwa and Kanematsu, 2018], enabling their practical application in everyday scenarios.The most relevant one was the use of an analytical amorphous track model of particle energy deposition at the nanometer and micrometer scale [Kiefer andStraaten, 1986, Chatterjee andSchaefer, 1976] which speed up the computation of microdosimetric quantities used as input for the MKM RBE models [Kase et al., 2007, Inaniwa andKanematsu, 2018].
In this work, we present a novel TOPAS Monte Carlo (MC) microdosimetric extension: MicrOdosimetry-based modelliNg for RBE ASsessment (MONAS).MONAS combines full MC simulations of microdosimetric spectra with clinically relevant microdosimetry-based radiobiological models for cell survival and dose-dependent RBE assessment.Furthermore, by utilizing the new scorers of MONAS, it is now possible to generate fully Monte Carlobased Look-Up Tables (LUTs) of radiobiological parameters α and β for RBE-based treatment plan optimization in clinical proton therapy.MONAS is based on the existing TOPAS microdosimetric extension [Zhu et al., 2019] which models the main microdosimetry detectors used in literature and scores the lineal energy distributions.The lineal energy y is defined as the energy deposited over the target volume mean chord length l, i.e. y = /l.This is the quantity of reference measured in experimental microdosimetry.Starting from the simulated y distributions, we implemented a novel scorer based on specific energy z, defined as the energy imparted over the mass m of sensitive volume, i.e. z = /m.The extension allows the user to calculate specific energy distributions at different microscopic scales.These z distributions are the building blocks for the microdosimetry-based radiobiological models implemented in MONAS.We included the GSM 2 and three clinically relevant MKM formulations: saturation corrected MKM (MKM-z*) [Kase et al., 2006], Double Stochastic MKM (DSMKM) [Sato and Furusawa, 2012] and the modified Stochastic MKM (mSMKM) [Inaniwa and Kanematsu, 2018].Therefore, MONAS predicts dose-dependent cell survival fraction and RBE specifically for the simulated radiation field.
MONAS simulates experimental microdosimetric distributions acquired with three different detectors, and from validated z-distributions, it calculates cell survival curves which can be directly compared with experimental in vitro data.Therefore, MONAS is the first full MC toolkit that allows the user to benchmark the physical input and the biological output of the radiobiological models with experiments.To show that, we compared the MONAS cell survival curves with experimental data from the Particle Irradiation Data Ensemble (PIDE) [Friedrich et al., 2013].We also calculated cell survival and RBE depth curves using all radiobiological models available in MONAS for a proton Spread-out Bragg-peak (SOBP), whose microdosimetric spectra were previously measured by us [Missiaggia et al., 2023] and used to validate TOPAS.
MONAS is also an accurate and fast tool to predict RBE in proton therapy treatment plans.In particular, in this study, we used the MONAS extension to generate full Monte Carlo-based Look-Up-Tables (LUTs) of radiobiological parameters α and β from monoenergetic proton beams.By combining MONAS LUTs and TOPAS MC's precision in tracking therapeutic protons within the patient's anatomy, MONAS allows the calculation of the RBE spatial distribution in a real patient's geometry.In particular, we determined the mixed-filed α mix and √ β mix , defined as the dose-averaged of single particle α and √ β extrapolated from MONAS LUTs, in each voxel of the scoring mesh.Consequently, we could predict the biological effectiveness of a real proton beam.As an example, we recalculated with TOPAS MC a real head and neck proton treatment, optimized with the Eclipse planning system (Varian Medical Systems, Palo Alto, CA) and delivered at the Dwoskin Proton therapy center (University of Miami).Our study demonstrated that by fully integrating the MONAS toolkit within the TOPAS MC code, we can accurately simulate microdosimetric spectra and use them as input for the microdosimetry base models, enabling us to estimate the effect of radiation on biological tissue in clinical conditions.The present work represents the first crucial step toward bridging the gap between microdosimetry and clinical applications.By providing a comprehensive microdosimetric analysis of the radiation field, we show how real-life clinical and experimental scenarios can greatly benefit, allowing clinicians and researchers to estimate the biological effect of radiation accurately.
The main contributions of the present paper are: (i) to introduce a Monte Carlo-based microdosimetric toolkit that allows estimating in vitro cell survival data and RBE distribution in realistic radiotherapy treatment plans with a full microdosimetric description of radiation.
(ii) to introduce a new methodology for utilizing the microdosimetry-based radiobiological models to evaluate RBE in real proton therapy treatments.

Microdosimetric distributions and their moments
Microdosimetry considers two quantities of interest: the specific energy z and lineal energy y [Zaider et al., 1996].
The specific energy z is the ratio between the energy imparted by ionizing radiation and the mass m of the sensitive volume, z = m .
The lineal energy y is the ratio between by ionizing radiation and the mean chord length of the sensitive volume l, y = l .
It is possible to relate the lineal to specific energies via the following relation z = l m y.By assuming a spherical site of density ρ t = 1 gcm −3 and radius r the relations between the y and z is: where 0.16 is the Gy − keV conversion coefficient.
The lineal energy y is the reference quantity in experimental microdosimetry, whereas the specific energy z is the main quantity of reference in microdosimetry-based radiobiological models.A key difference between lineal energy and specific energy is that the quantity refers to the energy imparted in a single event for y, or the energy imparted in any number of events for z.This implies that when considering the distribution of z, the ionization of more than one event must be considered.For this reason, most of the MKM and the GSM 2 models are based on the so-called multi-event specific energy distributions.
The multi-event distribution is defined starting from the n-event distribution f n (z) which defines the specific energy frequency distribution when exactly n events occurred.The nevent distribution can be computed as the n-fold convolution of the single-event specific energy distribution f 1 (z): As standard in microdosimetry, the first two moments of the single-event distribution play a crucial role and they are defined as The multi-event distribution is defined as [Bellinzona et al., 2021b, Zaider et al., 1996], where, by statistical independence of the n events they are assumed to follow the Poisson distribution with mean value λ n , so that the term multiplying f n in equation ( 6) is the probability of registering exactly n events.With further computations, [Zaider et al., 1996], it is possible to show that the first and second moments of the multi-event distribution are: Typically, z in equation 7 is identified with the absorbed dose D, yielding that the mean value of Poisson distribution of equation 6 is given by λ n = D z F .Further computation, [Zaider et al., 1996, Bellinzona et al., 2021b], shows 2.2.Microdosimetry-based RBE models 2.2.1.The Microdosimetric Kinetic Model and its generalizations The MKM is a radiobiological model that utilizes a system of differential equations to predict the survival fraction of irradiated cells.These equations describe the time evolution of the average number of DNA damages, which are divided into two types in the MKM: sub-lethal and lethal lesions.Sublethal lesions are DNA damages, typically DNA Double-Strand Breaks (DSB), that can be repaired by the cell, while lethal lesions are unrepairable and result in cell inactivation.To better align with biological data, the MKM postulates that the nucleus is partitioned into sub-units called domains so that the number of lethal and sublethal lesions is evaluated for each domain separately.Therefore, the probability of cell survival is estimated by taking into account all the domains into which the nucleus has been divided.
Since its original formulation in [Hawkins, 1994], the MKM has been widely generalized to include several endpoints and stochastic inter-and intra-cellular effects, [Hawkins, 1996, Hawkins, 2003, Inaniwa et al., 2010, Manganaro et al., 2017, Inaniwa and Kanematsu, 2018, Bellinzona et al., 2021b, Attili et al., 2022].We will refrain from providing an exhaustive explanation of all the MKM derivations.Instead, we will focus solely on recalling the formulas for cell survival fraction derived by the models.
Assuming that the reference radiation exhibits an LQ dependence as in equation ( 1) between the logarithm of the cell survival fraction and the imparted dose, with given parameters α X and β X , in its original formulation [Hawkins, 1994], the MKM predicts an LQ cell survival fraction of the form where β = β X coincides with the β of the reference radiation and α is where α 0 is the limit of the α parameter as the LET becomes increasingly small and zd,D is the dose-average specific energy defined in equation ( 5) computed over the cell nucleus domain.
A proposed correction to the MKM model, named saturation corrected MKM [Kase et al., 2006], aims to improve its alignment with heavy ion data, which exhibits the so-called overkill effect.This effect consists in a decrease of the RBE versus LET, for LET beyond approximately 150 keV/µm, following its initial raise for increasing LET [Kase et al., 2006].In particular, the dose-average specific energy zd,D is substituted with its saturation corrected value, [Zaider et al., 1996] where z 0 is the saturation correction parameter for modeling the over-killing effect due to high-LET radiation.
According to the saturation corrected MKM, the α parameters thus becomes whereas the LQ behavior (9) remains valid.
Although the saturated corrected MKM shows a better match with experimental data, it still does not include energy deposition variations both at the cell nucleus and at the cell population level.To include these effects, the DSMKM, [Sato and Furusawa, 2012], has been proposed.The cell survival fraction at nucleus scale S n as a function of nucleus-specific energy z n can be calculated from domain multi-event specific energy distributions as where the domain saturation-corrected specific energy is defined: . Considering the stochastic nature of z n , the survival fraction of cells irradiated by absorbed dose D, S(D), can be estimated by: The computational requirements for evaluating multi-event functions such as f (z d , z n /z d,F ) and f (z n , D/z n,F ) are non-negligible, particularly when combined with macroscopic beamtransport simulations that use Monte Carlo methods to model the irradiation conditions.Therefore, in [Sato and Furusawa, 2012] S n (z n ) formulation has been simplified by omitting the calculation of f (z d , z n /z d,F ), assuming instead that the saturation effect triggered by multiple radiation events to a domain is negligibly small, yielding where the fluence-and dose-averaged specific energy (z d,F , zd,D ) are calculated according to the equations ( 4), ( 5), respectively.The fluence-averaged saturation-corrected specific energy is defined z * , while the fluence-averaged saturation-corrected specific energy corresponds to equation ( 11).This formulation was called stochastic-MKM (SMKM).
The SMKM has been further simplified in [Inaniwa and Kanematsu, 2018] to speed up the computational time and to be implemented in treatment planning systems.The modified version of SMKM (mSMKM) is based on the assumption that, in charged-particle therapy, the domain-specific energy z d is in general delivered by a large number of low-energy deposition events, and the events inducing the saturation of complex DNA damages, i.e. events with z d > z 0 , are rare.Also, it is assumed that the specific energy imparted z n is sufficiently close to the macroscopic dose D. According to the mSMKM the cell survival fraction is calculated as follows: 2.2.2.The Generalized Microdosimetric Model (GSM 2 ) Recently, starting from the building assumptions of the MKM, a novel microdosimetry-based radiobiological model, GSM 2 , was presented, [Cordoni et al., 2021, Cordoni et al., 2022b, Cordoni et al., 2022a].GSM 2 aims at providing a fully probabilistic model which takes into account the effects of stochasticity in different aspects of radiation-induced damage, e.g., in the initial damage distribution as well as damage evolution.Similar to the MKM formulations, GSM 2 describes the time evolution of sub-lethal lesion x and lethal lesions y in a cell nucleus, which is divided into smaller sub-domains.Differently from the MKM, GSM 2 is able to describe the time evolution of the whole probability distribution of lesions rather than simple average values.
The survival fraction of a single cell due to the absorption of specific energy z n in the cell nucleus S n (z n ) is computed by GSM 2 as, [Cordoni et al., 2022b], where p X 0 (x 0 |z n ), p Y 0 (0|z n ) represent the probability that the domain suffers x 0 sub-lesions and 0 lethal lesions, respectively, N d is the number of domains in which the cell-nucleus is divided and C(x 0 ) are weighting terms representing the probability that x 0 sub-lethal lesions are repaired.The terms p X 0 (x 0 |z n ), p Y 0 (0|z n ) and C(x 0 ) can be calculated as follows: with p X (x|kz) and p Y (y|λz) being Poisson distribution of average kz and λz, respectively.Notable enough, it has been shown in [Cordoni et al., 2022b, Missiaggia et al., 2022], that the distributions p X 0 (x|z n ) and p X 0 (y|z n ) can deviate from a Poisson distribution, with particular emphasis on clear differences for sufficiently high LET and doses.Further, r, a, b, k, λ are model parameters.In particular, r is the repair rate, a is the rate at which a sub-lethal lesion becomes a lethal lesion and b is the pairwise interaction rate of two sublethal lesions to become a lethal one.At last, λ and κ are the average yields of lethal and sublethal lesions, respectively, per unit Gray.
Also, another remarkable property of the predicted cell survival fraction ( 18), is that it exhibits an LQ behavior at low doses, naturally straightening to a purely linear trend at high doses, [Cordoni et al., 2022b, Missiaggia et al., 2022].This pattern matches empirical evidence collected in cell survival experiments.
At last, as for MKM formulations, the cell survival fraction S(D) of a cell population due to the absorbed dose D can be estimated according to the equation ( 14).

TOPAS
The microdosimetry-based radiobiological extension presented in this work extends the TOPAS MC toolkit [Perl et al., 2012].TOPAS is an easy-to-use interface to the Geant4 Simulation toolkit [Agostinelli et al., 2003] allowing both medical physicists and researchers to make Monte Carlo simulations without the necessity of advanced coding knowledge.In [Zhu et al., 2019] the microdosimetric extension of TOPAS has been implemented, allowing to score of lineal energy y with three different types of detectors: (i) spherical Tissue Equivalent Proportional Counter (TEPC), (ii) a cylindrical TEPC (also known a mini-TEPC) and (iii) Silicon on Insulator (SOI) microdosimeter.These detectors are the reference detectors for microdosimetry, [Missiaggia et al., 2020, Missiaggia et al., 2023, De Nardo et al., 2004, Bianchi et al., 2022, Bradley et al., 2001, Debrot et al., 2018].The microdosimetric extension offers to the user the possibility to save microdosimetric spectra (yf (y) and yd(y)) and the relative average quantities (y F and y D ) for each detector, including the contribution of the particle species of the radiation field.The lineal energy distributions obtained via the TOPAS microdosimetric extension have been benchmarked with experimental data, showing good agreement, [Zhu et al., 2019, Missiaggia et al., 2023].
2.3.1.MONAS MONAS extension starts from the original lineal energy scorer to provide a further toolkit that calculates specific energy z, as introduced in Section 2.1.Based on specific energy microdosimetric spectra, MONAS predicts cell survival fraction and RBE using the different MKM formulations as well as GSM 2 radiobiological model.The workflow of MONAS is depicted in Figure 1.
In addition to the parameters of the lineal energy scorer [Zhu et al., 2019], new optional parameters were implemented to activate the cell-survival and RBE calculations module according to MKM and GSM 2 : GetRBEWithMKModel, GetRBEWithGSM2.The user can further choose one or more MKM formulations available by setting the value of the string parameter MKMCalculation equal to i) "MKM-z* " for the saturation corrected MKM, ii) "mSMKM " for the modified-SMKM or iii) "DSMKM " for the double stochastic MKM.By default, the saturation corrected is set.The cell-survival S is computed as a function of macroscopic absorbed dose D, given as input parameter by the user setting the parameter SurvivalDoses; thus RBE values are calculated as a function of S(D) as follows where α X and β X are the Linear-Quadratic coefficients of photon reference radiation.Modelspecific radiobiological parameters can be set both for MKM and GSM 2 separately, including the reference radiation coefficients for the RBE calculations.Table 1 summarizes all MONAS parameters and their default values.Further details about the evaluation of radiobiological parameters will be given in Section 2.4.1.
The output files include the values set by the user both for the radiobiological model and reference radiation, specific energy spectra both for cell nucleus and domain, cell survival for the macroscopic doses specified, and RBE as a function of cell irradiation dose in ASCII format.The MONAS extension is an open-source code available on GitHub [Cartechini, 2023].

Specific Energy Spectra
Parallel to the default lineal energy scorer, MONAS includes the calculation of specific energy quantities converting single-event lineal energy y into singleevent specific energy z 1 according to the equation (2).Single-event and multi-event specific energy spectra are calculated and used for cell survival and RBE evaluation, as described in Section 2.2.By setting the boolean parameter SaveSpecificEnergySpectra, the user can save in ASCII format the single-event and multi-event distributions calculated on cell domain and cell  nucleus: , where the subscript x can be either d for the domain and n for the nucleus distributions.Since the n-fold convolution for the multi-event calculation is time-consuming, especially for high doses when the number of convolutions increases, a Monte Carlo approach for evaluating the multi-event distribution has been specifically implemented according to the following workflow: (i) the number of tracks k which deposit energy on sensitive volume is generated from a Poisson distribution with mean value λ n = z n /z d,F for cell domain and λ n = D/z n,F for nucleus, respectively.Then, (ii) k single-event specific energies z 1 are sampled from the single-event probability distribution f 1 (z) and summed up to obtain the total specific energy z tot = k i=1 z 1,i deposited in the target; (iii) the multi-event probability distribution f (z, λ n ) both for cell domain and the nucleus is thus constructed by iterating steps (i) and (ii) N-times according to the parameter SetMultieventStatistic.A scheme of the algorithmic construction described is depicted in figure 2.
Determine the number of hits, k, in the domain (nucleus) from a Poisson distribution of mean value z n /z F (D/z F ) • For each hit, the specific energy z i is obtained by sampling the single-event probability distribution f 1 (z) • The multi-hit specific energy is update z tot,j = z tot,j + z i for i = 0,…,k For j =1,…,N Accumulate the z tot,j on a frequency histogram Normalize the frequency histogram to obtain the multi-event probability distribution f(z,l n ) for cell domain (nucleus) Figure 2: Flowchart of multi-event probability distribution calculation with a Monte Carlo approach by sampling distributions via random number generators.Steps inside curly brackets are the body of the two for loops iterating on the number of multi-event statistics and cell domain (nucleus) particle hits.

Cell-survival and RBE
The main feature of the MONAS extension is the prediction of cell survival curves and dose-dependent RBE after irradiation with ion beams.To prove the accuracy of the toolkit, we determined the radiobiological parameters specific to each model by comparing the survival predictions with experimental in vitro data on the Human Salivary Gland (HSG) cell line.Once the model parameters were determined, we used MONAS to predict RBE distribution in two relevant cases for proton therapy: passively scattering proton Spread-Out Bragg Peak (SOBP) [Missiaggia et al., 2023, Tommasino et al., 2019] and head and neck proton therapy irradiation delivered at the Dwoskin Proton therapy center (University of Miami) and optimized with the Eclipse treatment planning system (Varian Medical Systems, Palo Alto, CA).

Radiobiological parameters
The radiobiological models implemented in this work are based on a variable number of free parameters which are independent of the radiation type, but they are only cell-line dependent.
These parameters were estimated by fitting the models with in vitro cell-survival experimental curves available in the literature and measured for a specific irradiation condition.It is worth stressing that, previous works [Kase et al., 2006, Sato and Furusawa, 2012, Inaniwa and Kanematsu, 2018] reported MKM parameters for HSG cell line.Nonetheless, since the physical estimation of the radiation field is different from what was implemented in MONAS, the biological parameters have been fitted to ensure the highest reproducibility of in vitro cell survival experiments.In the present work, we characterized the radiation field in terms of lineal energy and specific energy spectra by exploiting the TOPAS toolkit, a condensed history Monte Carlo algorithm.On the contrary, a different methodology was used both for the DSMKM and for the modified-SMKM [Sato andFurusawa, 2012, Inaniwa et al., 2013].In fact, DSMKM exploits a combination of microdosimetric Monte Carlo simulations with PHITS code [Sato et al., 2013] in 1 µm volume for the cell domain [Sato et al., 2009, Sato et al., 2006], while a Fermi function of macroscopic LET is used for evaluating the specific energy probability distribution on cell nucleus [Sato and Furusawa, 2012].The modified-SMKM, instead, exploits analytical amorphous-track models to describe the radial dose distribution of the ion track and, thus, to compute the specific energy deposited on the cell domain and nucleus [Chatterjee and Schaefer, 1976, Kiefer and Straaten, 1986, Kase et al., 2006, Inaniwa et al., 2010].Due to the intrinsic difference in the description of radiation energy deposition in the sensitive volume between TOPAS and the previous work, we re-calculated the model parameters.A systematic comparison between Monte Carlo condensed history (e.g.TOPAS MC) and track structure algorithm for the calculation of specific energy spectra is out of the scope of this work.
Therefore, the parameters for the MKM and GSM 2 were determined to reproduce in vitro experimental data of HSG cells.Experimental cell survival curves were taken from the Particle Irradiation Data Ensemble (PIDE) [Friedrich et al., 2013] in which a large amount of cell survival data are systematically collected and analyzed as a function of particle LET, cell line, and reference radiation.Regarding the estimation of specific energy spectra, we simulated with TOPAS MC (v3.7) a water sphere placed in a vacuum world to avoid particle energy loss outside the sensitive volume.The sphere was irradiated by a mono-energetic 3 He (10.2 MeV/u and 4.89 MeV/u) and 12 C (126 MeV/u and 126 MeV/u) ion beams as reported in the experiments [Furusawa et al., 2000, Friedrich et al., 2013].The beam was modeled as the Environment particle source available in TOPAS.It creates an isotropic, uniform radiation field enclosing the water sphere.The default TOPAS physics list was used [Jarlskog and Paganetti, 2008].Then, we scored the specific energy spectra in the sensitive volume.
2.4.2.Cell survival and RBE for a proton Spread-out Bragg-peak Recently in [Missiaggia et al., 2023] a systematic characterization of the radiation field produced by passively scattered proton SOBP generated by a pencil beam of 148 MeV has been presented.Microdosimetric spectra were acquired with spherical TEPC both in-field and out-of-field.The same work also shows a good agreement between TOPAS microdosimetric simulations and experiments.
Starting from validated lineal energy spectra, we calculated cell survival, and RBE along the beam axis using the MONAS code.As described in [Missiaggia et al., 2023], we simulated the Far West Technologies LET-1/2 spherical TEPC with an active volume made of pure propane gas (C 3 H 8 ) at an operative pressure such that it is equivalent to a tissue sphere of 2 µm diameter.The default physics list was used and the particle production cut was set according to the microdosimetric extension [Zhu et al., 2019].About 10 7 primaries were simulated for each position along the beam axis.200 kV X-ray was chosen as reference radiation for the HSG cell line, with α X = 0.19 Gy −1 and β X = 0.05 Gy −2 [Kase et al., 2006].
2.4.3.Cell survival and RBE assessment in a patient case Predicting the cell survival fraction and RBE within the patient's anatomy involved a two-step process.First, we utilized the MONAS extension to create lookup tables (LUTs) of radiobiological parameters α and β for each model and proton beam energy.Then, we incorporated the MONAS LUTs into the TOPAS Monte Carlo particle transport algorithm to score the survival fraction and RBE in each voxel of the patient.To construct the LUTs, we irradiated a sphere of 1 µm of radius made of water with monoenergetic proton beams at different energies (from 0.1 MeV to 300 MeV) and we scored the cell survival fraction with the MONAS extension.We fitted each survival curve with the LQ model and we determined the α and β coefficients.Using the TOPAS particle transport algorithm, we created a novel scorer to calculate the mixed filed α mix and β mix for each voxel of the patient scoring mesh starting from the LUTs of monoenergetic beams.In particular, for each step s of primary and secondary protons inside the voxel v, we registered the kinetic energy of the particle (E v,s ) and the energy deposited along the step plus the energy released to secondary δ-rays (dE v,s ) [Cortés-Giraldo and Carabe, 2015].As standard, the α mix and β mix are thus calculated as a weighted sum of α and β for the monoenergetic beam, [Zaider and Rossi, 1980], To show the applicability of this approach in real patient irradiation, we simulated in TOPAS a head and neck treatment plan, optimized with the Eclipse treatment planning system (Varian Medical Systems, Palo Alto, CA) at the Dwoskin Proton therapy center (University of Miami).The plan is composed by 2 coplanar fields at 30 and 60-degree gantry angles and a third noncoplanar at 30 degrees gantry and couch angle at 300 degrees.All fields employed a range shifter of 57 mm of water equivalent thickness.A uniform biological dose of 60 Gy(RBE) in 30 fractions was prescribed to the target using a constant RBE equal to 1.1.
We estimated the cell survival fraction and the dose-dependent RBE with the approach described above.The radiobiological parameters for the reference radiation were α X = 0.19 Gy −1 and β X = 0.05 Gy −2 for HSG cell line.
To compare our findings, we repeated the same analysis using the mSMKM based on the amorphous track description of radiation energy description at the microscopic scale [Inaniwa and Kanematsu, 2018], which is the version currently used in carbon ion therapy.To generate the dose-averaged specific energies per event zd,D and z * d,D imparted to the domain, and to the cell nucleus, zn,D , we employed the Survival toolkit code [Manganaro et al., 2018].The resulting cell survival fraction and RBE are as described in [Inaniwa and Kanematsu, 2018]. .Experimental data from [Furusawa et al., 2000], tabulated in the PIDE.

Radiobiological parameters
Table 2 reports the MKM and GSM 2 parameters that give the best fit for the in vitro cell-survival data of the HSG cell line when irradiated with 12 C and 3 He ion beams.Figure 3 shows the experimental cell survival fraction, as taken from the PIDE dataset, [Friedrich et al., 2013], for the HSG cell line, compared to the corresponding predicted cell survival curves using the four radiobiological models with parameters as given in Table 2.All radiobiological models are in good agreement with the HSG in vitro cell survival curves for both carbon and helium ions.
Regarding the fitted MKM parameters, given a new description of the physics of the radiation field employed in this study, all parameters have been recalibrated due to a discrepancy between cell survival experimental results and prediction using parameter values reported in the original papers.In particular, in [Kase et al., 2006], the MKM-z* parameters extrapolated directly from microdosimetric measurements of y D and from in-vitro HSG cell survival data are reported to be α 0 = 0.13 Gy −1 , β = 0.05 Gy −2 .Absolute values of α 0 and β do not agree with our fit, nevertheless, the α/β ratio and domain radius are consistent with [Kase et al., 2006].Further, DSMKM and mSMKM parameters are different from the originals, [Sato andFurusawa, 2012, Inaniwa andKanematsu, 2018], but nevertheless, the difference in absolute value is moderate.
At last, concerning GSM 2 , coherently with [Missiaggia et al., 2022], only a, b, and r parameters were fitted whereas domain and cell nucleus radius were set a priori.This is done to avoid overfitting.
Table 2: MKM and GSM 2 model parameters, as reported in Section 2.2, for Human Salivary Gland cell line.Experimental data are taken from the PIDE dataset, [Friedrich et al., 2013].

Specific energy spectra
The novel specific energy scorer implemented in this work allows calculating the single-and multi-event probability distributions in equation ( 6), both at the domain and cell nucleus scales, that is, in the order of 1 micron and 10 microns, respectively.Figure 4 shows the single-event zf 1 (z) and multi-event zf (z, z n ) distributions computed at different averages doses z n =0.1, 1 and 30 Gy for a cell domain radius equal to 0.46 µm, chosen as the representative domain sizes for the DSMKM and the mSMKM.These spectra were calculated by locating the TEPC active volume at the center of the SOBP proton beam described in Section 2.4.2.
The single-event distribution zf 1 (z) preserves a similar shape as the corresponding yf (y) distribution, as could be expected from equation (2).On the contrary, the shape of the multievent distribution changes significantly as a function of the dose delivered to the cell domain z n .In particular, to lower doses, e.g.z n =0.1 and 1 Gy, it corresponds a higher probability of null energy deposition in the domain due to the case of no tracks hitting the target.This phenomenon is evident by the high peak at z = 0 Gy.As the dose z n increases, since the average value of Poisson distribution is proportional to z n , the probability to score zero tracks vanishes, and consequently the peak at null z disappears.All three multi-event distributions have an average value equal to z n , and furthermore, at higher z n the distribution is uni-modal and peaked at z = z n with a Gaussian-like shape.
The zf (z, D) distribution at the cell nucleus scale, with a radius of 8 µm, as a function of macroscopic dose D =1, 5, 15 Gy is shown in Figure 5.As before, the spectrum is generated by an SOBP proton beam.The multi-event distribution has a uni-modal Gaussian-like shape peaked at z = D even at a low dose of D = 1 Gy.

Cell survival and RBE for a proton Spread-out Bragg-peak
Cell survival fraction and RBE were calculated using MKM formulations and GSM 2 along the beam axis using experimentally validated microdosimetric spectra, [Missiaggia et al., 2023], as described in Section 2.4.2.The physical 3D dose distribution was simulated in the water phantom using a voxel size of 1 × 1 × 1 mm 3 and normalized to 1.8 Gy at the center of SOBP to obtain a biological dose of 2 Gy(RBE) for constant RBE equal to 1.1.
Figure 6 panel (a) shows the depth survival curve compared to the physical dose.All models predict a similar trend: the survival is higher at the entrance, around 0.7, and drops to a minimum in the SOBP, where the survival remains constant at around 0.5.At the distal edge of the SOBP, the survival fraction raises again to a similar value as the entrance channel.However, quantitative differences between the considered models are appreciable.The MKM-z* and mSMKM formulations predict cell survival fractions systematically lower than GSM 2 in all regions.The DSMKM shows comparable survival fraction values with GSM 2 in the entrance up to the mid-SOBP, while in the distal region, GSM 2 predicts a slightly higher cell survival.Dose-dependent cell survivals were further used to calculate RBE according to equation ( 22).
All the models predict an RBE higher than 1 for all depths, Figure 6 panel (b).In addition, all models show a constant RBE trend as a function of penetration depth with a value between 1.1 and 1.2 in the entrance channel and a sharp increase in the distal edge of the field.Table 2 reports the survival fraction and RBE values, along with statistical errors, at three water depths: entrance (34 mm), middle SOBP (116 mm), and distal edge (136 mm).0.77 ± 0.07 1.47 ± 0.08 MKM-z* 0.70 ± 0.04 1.86 ± 0.07 mSMKM 0.71 ± 0.04 1.82 ± 0.07 DSMKM 0.72 ± 0.07 1.81 ± 0.09 Table 3: Cell survival fraction, S, and RBE calculated with MKM formulations and GSM 2 for HSG cell line.The values are reported at three relevant water depths of the proton SOBP: entrance (34 mm), mid-SOBP (116 mm), and distal penumbra (136 mm).

MONAS application to a patient case
In order to estimate a 3D spatial distribution of cell survival and RBE for a patient case, we calculated lookup tables (LUTs) of radiobiological parameters α and β for protons specific to the HSG cell line.
We incorporated the LUTs within the TOPAS Monte Carlo algorithm and simulated a treatment plan, calculating the 3D spatial distributions of dose-dependent cell survival and RBE.To benchmark our results, we also included the mSMKM with amorphous track (mSMKM-AT), a validated radiobiological model for Carbon ion therapy, [Inaniwa and Kanematsu, 2018].The 2D spatial distribution of cell survival fraction is plotted in Figure 7 for GSM 2 (panel (c)), and mSMKM-AT (panel (a)) as a benchmark to the current state-ofthe-art of MKM clinical application (tuned for carbon ion therapy).The two models agree in the shape of S(D) distribution where a minimum (color map hot region) of the cell survival fraction is estimated in the Clinical Target Volume (CTV) and a few millimeters in the surrounding volume.By comparing the monodimensional depth-survival curves (panel (e)), large differences are evident among the models.MONAS-based MKM formulations predict the lowest S(D) values (0.47 ± 0.03 on average in CTV), while GSM 2 and mSMKM-AT predict an average cell survival fraction of 0.51± 0.07 and 0.55 ± 0.04, respectively.
The RBE 2D map displayed in panels (b) and (d) exhibits similar patterns: both the mSMKM-AT and GSM 2 are characterized by a notable increase in RBE values in the distal region of the therapeutic field, specifically at beam angles of 30 and 60 degrees.mSMKM-AT 2D distribution shows a wider dark red region (RBE values above 1.5) compared to GSM 2 in the distal part of the CTV, while in the patient's entrance mSMKM-AT map has lower RBE values (blue colormap).As shown in the depth-RBE curve (panel (f)), MKM models from MONAS code predict RBE values systematically higher than GSM 2 and mSMKM-AT both at the beam entrance and in the distal region (from 1.3 to 1.6).GSM 2 model predicts RBE values lower than the MKM formulations from around 1.1 at the entrance, up to 1.5 downstream of the CTV.The same holds for mSMKM-AT which shows RBE values around 1 at the entrance, exceeding GSM 2 in distal reaching RBEs of around 1.7.
Although the inter-model differences in absolute values of the survival S and RBE, both MONAS-based and mSMKM-AT models show the same spatial distribution characteristics.Inside the target volume S and RBEs remain constant and sharply increase at the distal edges of the fields where most of the organs at risk are located.

Discussion
In this work, we presented the MONAS toolkit, a TOPAS MC extension that combines accurate Monte Carlo simulations of microdosimetric distributions with microdosimetry-based models for predicting cell survival and RBE.MONAS provides a powerful tool to bridge the gap between a microdosimetric description of a radiation field, the biological evaluation of radiation effects, and its implementation in the clinical practice of proton beam therapy.
To show MONAS applications, we presented two examples: i) the radiobiological characterization of a SOBP, and ii) the calculation of the RBE spatial distribution for a real patient's treatment with protons.In Figure 6, we showed that MONAS can estimate the cell survival fraction and RBEs from experimentally validated spectra, predicting a trend consistent with previous works [Kase et al., 2006, Tran et al., 2017].In this context, it is important to emphasize that most of the research paper on protons and heavier ions focuses on RBE 10 , which is the ratio of reference radiation and ion dose giving the 10% survival fraction [Bianchi et al., 2020, Conte et al., 2020, Debrot et al., 2018, Tran et al., 2017].The RBE 10 represents a universally accepted radiobiological endpoint, allowing thus a comparison of different radiobiological experiments.However, when considering clinical applications, it becomes more meaningful to examine a dose-dependent RBE, i.e. the RBE associated with the dose specifically absorbed by each voxel of the treatment plan.This critical information enables direct estimation of the treatment plan's biological effectiveness.The presented MONAS extension, which is directly linked to TOPAS, offers a rapid assessment of the MC-based dosedependent RBE, which can have an impact both in clinical scenarios and in experimental radiobiology.
Using the complete range of microdosimetric spectrum as the input for radiobiological models emphasized the need to recalibrate the model parameters, as outlined in Table 1, in order to accurately fit the MONAS cell survival curves with the existing experimental data.Figure 3 shows a good agreement between MONAS predictions and the experimental in vitro data at different irradiation conditions.Furthermore, the new set of parameters is consistent with parameters already published for the MKM [Kase et al., 2006, Sato and Furusawa, 2012, Inaniwa and Kanematsu, 2018] in terms of the order of magnitude of the cell domain and nucleus radii, as well as the α-β absolute values.The MONAS toolkit thus further allows testing the robustness of the model parameter against experiments.In MONAS, all microdosimetric values are evaluated using the entire microdosimetric distribution (Figures 4-5), which can be validated by direct comparison with the experimental y-spectra measured with commercial microdosimetric detectors.As a result, the MONAS toolkit provides a comprehensive framework for validating the entire radiobiological workflow, encompassing the simulation of the radiation field's physical properties and the direct estimation of its biological impact.This is achieved by utilizing microdosimetric quantities and cell survival experiments that are directly calibrated against experimental data, ensuring accuracy and reliability in assessing biological effectiveness.
The second main application of MONAS is the RBE evaluation in a real patient case.In Figure 7, we show the 3D spatial distribution of cell survival fraction and RBE values calculated for proton treatment of a head and neck tumor.The results for a complex beam and patient's geometry follow the same trend exhibited in Figure 6 from the SOBP: RBE hot spots are observed at the distal edge of the field, where most of the organs at risk are located (e.g.salivary glands, oral cavity, optical nerves).We used the MONAS extension to generate LUTs of α and β coefficients for monoenergetic proton beams and eventually asses the biological effectiveness of a mixed-energy beam averaging α and β values.The RBE evaluated in the present paper is calculated considering only the contribution from primary and secondary protons, but it can be extended to include other secondary ions by adding LUTs for alpha, carbons, and oxygen ions.This could increase the accuracy of the estimated RBE values, especially at the field-edge and out-of-field, where a non-negligible contribution to the radiation field is given by high-LET radiation such as alpha particles [Missiaggia et al., 2023].The application of MONAS in a real patient case represents a significant advancement in the integration of a detailed microdosimetric description of the radiation field into treatment planning systems.These LUTs can be incorporated into a particle tracking Monte Carlo toolkit to estimate a three-dimensional map of cell survival fraction and RBE within a mixed-field radiation treatment plan.
While the use of LUTs does not fully address the computational challenge of simulating the complete microdosimetric distribution in all voxels, we propose a fundamental step into the usage of microdosimetric information in proton and ion therapy.LUTs have been generated by simulating the entire microdosimetry distribution of monoenergetic beams, which were then used as input for the radiobiological models to estimate the cell survival fraction.Therefore, the intrinsic stochastic nature of energy deposition at the micron scale is taken into account in the final LUT values.The use of pre-calculated α and β values that are loaded during the particle simulation in patients does not increase the treatment plan computation time.In fact, the time required by TOPAS MC to calculate α mix and √ β mix in each voxel of the patient is comparable to a similar volumetric scorer already implemented in TOPAS (e.g.ProtonLET [Cortés-Giraldo andCarabe, 2015, Granville andSawakuchi, 2015]).This feature is of great significance in the proposed methodology because it enables the inclusion of our LUTs in clinical Monte Carlo treatment planning systems without causing a slowdown in computational time for treatment plan calculations.
Our findings reveal that there is a clear inter-model variability in the absolute values of cell survival and, thus in RBE values observed both in SOBP characterization and in the patient simulation, as seen in Figures 6-7.Figures 7 panels (e)-(f) exhibit significant differences between the MKM models and GSM 2 in the CTV region, being further the mSMKM-AT closer to the GSM 2 prediction.Differences in the models stem from the fact that they all have different foundations and include different stochasticities.In addition, deviations in the predictions of RBE have already been verified in the literature when comparing the MKM with the LEM, [Monini et al., 2019, Bellinzona et al., 2021b].Furthermore, the predictions of the implemented models show similar patterns, with GSM 2 predicting a higher survival probability in the tumor compared to all MKM versions.The aim of this study is to provide a proof-of-principle of MONAS key features without making any assumptions about the model's parameters.The strength of our toolkit lies in the ability to customize each model's parameters and cell linespecific settings.However, we provided a set of potential parameters that yield radiobiological results consistent with the in vitro data presented in Figure 3.Further investigations will be conducted to calculate the model parameters that best fit a larger sample of experimental cell survival curves at various irradiation conditions.This analysis highlights the importance of determining precise model parameters that have a significant impact on the absolute values of cell survival fractions and, consequently, RBE.Furthermore, all models predict an RBE that significantly differs from the constant value used in the clinic.GSM 2 and the DSMKM calculate an RBE value close to 1.1 only in the entrance channel, while all the models show a sharp increase in the RBE in the distal region.This increase in the RBE could have an impact on the organs at risk located right after the tumor, for which a significant underestimation of RBE could lead to toxicities.
Microdosimetry is proving to be extremely valuable in clinical applications for several compelling reasons, surpassing the traditional use of LET values and providing a more reliable and experimentally measurable, as demonstrated in numerous campaigns conducted over the years, [Kase et al., 2006, Tran et al., 2018, Missiaggia et al., 2020, Bianchi et al., 2020, Lee et al., 2021, Missiaggia et al., 2023, Magrin et al., 2023], radiation quality description.Moreover, microdosimetry naturally includes the geometry of the sensitive volume in the measured spectra, providing information that not only considers the stochastic nature of energy deposition but also incorporates specific geometric considerations of the studied volume.This eliminates any potential uncertainty in interpreting the analyzed spectra and quantities.LET often fails to specify whether the value represents the track average or dose average.Additionally, different LET scorers are commonly used, and including only primary particles or both primary and secondary particles in the scorer can completely change the physical significance of the resulting values and the estimation of the biological effect [Grassberger and Paganetti, 2011, Bellinzona et al., 2021a, Kalholm et al., 2021].This becomes particularly problematic in the out-of-field regions, where the radiation field is mixed and short-track particles can play a significant role, and LET is known to be a bad descriptor of the radiation quality, [Grün et al., 2019].On the other hand, microdosimetry offers comprehensive information that allows for a precise assessment of all the physical processes involved, enabling a more accurate estimation of the biological effects.
Therefore, radiotherapeutic clinical practice can greatly benefit from consistently and rigorously evaluating treatment plans based on microdosimetry.In the past, the computational effort required for calculating microdosimetry quantities in treatment plans has been a hindrance to its practical implementation in the clinic.However, with advancements in computational power and the effective simplification methods demonstrated in this study, clinicians can now improve the prescribed treatment plans by incorporating microdosimetric considerations allowing for a more robust estimation of the plan's biological effectiveness.
MONAS wraps the already published TOPAS microdosimetric extensions to evaluate the single-and multi-event specific energy (z) distributions at different micrometric scales.Full microdosimetric distributions are then used as input for both MKM and GSM 2 models.This approach showed intrinsic differences in microdosimetric radiation characterization with respect to the amorphous track structure model used in the latest MKM formulations.Therefore, we recalculated the model parameter which best fit the radiobiological experiments for the HSG cell line.To show the main MONAS applications, we reproduced experimental microdosimetric spectra from a passively scattered SOBP.We used the MONAS code to assess cell survival fraction and RBE as a function of proton penetration depth.Our findings are consistent with the well-known RBE trend which presents a steep increase in the distal edge of the field.Furthermore, we were able to assess the high inter-model variability on the absolute RBE valuesthus quantifying a radiobiological uncertainty in proton plans in addition to other physical uncertainties.
The applicability of the MONAS toolkit can be further extended to a radiobiological analysis of treatment plans.We showed that it is possible to generate radiobiological parameter look-up tables which can be combined with the Monte Carlo toolkit for computing RBE maps on patients and therapeutic beam geometries.We showed an example of cell survival and RBE predictions on a real head and neck proton therapy plan delivered at the Dwoskin Proton therapy center at the University of Miami.The results have been compared to the mSMKM model based on the amorphous track structure model, as recently developed in [Inaniwa and Kanematsu, 2018], which represents one of the most used models in carbon ion therapy.Despite the variability in RBE absolute values, all the models showed a reasonable RBE trend as a function of beam penetration depth.Further investigation will focus on a full a biological robust optimization, based on the above mentioned uncertainty.
In conclusion, the MONAS extension offers a comprehensive microdosimetric framework for assessing the biological effect of radiation in both research and clinical environments.MONAS could be a key tool to include a detailed microdosimetric description of radiation field into treatment planning systems for variable RBE calculations.

Figure 1 :
Figure 1: MONAS workflow.Monoenergetic microdosimetric spectra are simulated using TOPAS MC to predict cell survival fraction and RBE with different radiobiological models.The predicted cell survival is thus used to construct α and β look-up tables.Then, voxel-based energy depositions are scored and used to reconstruct α and β parameters for mixed fields, from which RBE distribution in clinical patients is obtained.

Figure 4 :
Figure 4: Single-event (a) and multi-event (b) specific energy probability distributions calculated on cell domain of radius 0.46 µm.The spectra were scored by locating the sensitive volume at the center of a proton SOBP.

Figure 6 :
Figure 6: Radiobiological characterization of a proton SOBP.Cell survival fraction (panel (a)) and dose-dependent RBE (panel (b)), as a function of depth in water, predicted by DSMKM (black), GSM 2 (orange), MKM (light blue), and mSMKM (green).The SOBP region is depicted in a shaded grey.The horizontal dashed black line in panel (b) correspond to RBE 1.1, which is the standard value used in clinics.

Figure 7 :
Figure 7: Panels (a) and (b) show the 2D spatial distributions of cell survival and RBE calculated with mSMKM-AT on the patient case, respectively.The same patient slice and distributions are shown for the GSM2 model, representative of the MONAS prediction (panels (c) and (d)).The horizontal dotted white line indicates the position of the line curves plotted in panels (e) and (f) for cell survival and RBEs, respectively.The red vertical band in the latter plots shows the position of the target region.

Table 1 :
Summary of new input parameters for survival, RBE, and quality factor scorers.The parameter types are indicated according to the TOPAS syntax (b stands for boolean, sv for string vector, u for unitless double, i for integer).Default model-specific biological parameters refer to the HSG cell line.