www.nature.com/scientificreports OPEN Radiation dose estimation with time‑since‑exposure uncertainty using the γ ‑H2AX biomarker Dorota Młynarczyk 1*, Pedro Puig 1,2, Carmen Armero 3, Virgilio Gómez‑Rubio 4, Joan F. Barquinero 5 & Mònica Pujol‑Canadell 5 To predict the health effects of accidental or therapeutic radiation exposure, one must estimate the radiation dose that person received. A well-known ionising radiation biomarker, phosphorylated γ -H2AX protein, is used to evaluate cell damage and is thus suitable for the dose estimation process. In this paper, we present new Bayesian methods that, in contrast to approaches where estimation is carried out at predetermined post-irradiation times, allow for uncertainty regarding the time since radiation exposure and, as a result, produce more precise results. We also use the Laplace approximation method, which drastically cuts down on the time needed to get results. Real data are used to illustrate the methods, and analyses indicate that the models might be a practical choice for the γ -H2AX biomarker dose estimation process. Ionising radiation is currently used for a variety of purposes, such as industrial radiography, energy production, and health diagnosis and treatment. This results in an increased risk of radiological accidents. When a radiologi- cal accident occurs, it is important to assess the dose of radiation absorbed by people affected to help decision- making measures and whether it is necessary to help with medical care. Biodosimetry refers to the use of several biological markers to quantify the exposure to radiation of any person suspected of being exposed. The present study focuses on the γ -H2AX protein, which is an accepted biomarker of dose exposure to ionising r­ adiation1. When a cell is exposed to radiation, a wide variety of DNA damage can occur, including double-strand breaks (DSBs), which literally means that the two strands of the double helix break in close proximity. The cell has DNA damage response mechanisms, that they are activated in the face of this type of disruption. A critical modification is that in presence of DSBs the histone variant H2AX is phosphorylated (γ -H2AX) near the break site. This modification is extended to many megabases of the chromatin allowing to be detected microscopically in form of nuclear foci by immunostaining ­techniques2, as seen in Fig. 1. Those visible fluorescent foci can then be scored manually or automatically by a computer program. Because foci yield is proportional to the absorbed radiation dose, the results from several hundred cells can be used for dose estimation a­ nalysis1,3. The first part of the dose estimation method involves the calibration process. Samples of human peripheral blood lymphocytes are irradiated under controlled conditions at predetermined doses. Then the number of foci counts in peripheral blood mononucleated cells (PBMCs) is collected and the results are fitted to regression curves or surfaces, giving a calibration model. To determine the absorbed radiation dose by an exposed person, the number of γ -H2AX foci in the blood sample of that person must be evaluated. Then an inverse regression ­method4,5, using the calibration model and the foci counts in the patient’s sample, allows to estimate the dose received by the patient. It is known that the number of DSBs depends on the amount of radiation received by cells; the more radiation absorbed, the greater the biological damage. However, the γ -H2AX foci are influenced not only by dose but also by the time from ­exposure6. As studied before, the γ -H2AX level reaches the maximum 30 min after exposure and then gradually decreases over hours as the cells repair the damage, returning to normal levels within 24–48 ­h7,8. For this reason, in the case of a radiological accident the utility of γ -H2AX assay is restricted to rapid usage. 1 Departament de Matemàtiques, Universitat Autònoma de Barcelona, Bellaterra, Spain. 2Centre de Recerca Matemàtica, Bellaterra, Spain. 3Departament d’Estadística i Investigació Operativa, Universitat de València, València, Spain. 4Department of Mathematics, School of Industrial Engineering, Universidad de Castilla-La Mancha, Albacete, Spain. 5Departament de Biologia Animal, Biologia Vegetal i Ecologia, Universitat Autònoma de Barcelona, Bellaterra, Spain. *email: dorotaanna.mlynarczyk@uab.cat Scientific Reports | (2022) 12:19877 | https://doi.org/10.1038/s41598-022-24331-1 1 Vol.:(0123456789) www.nature.com/scientificreports/ Figure 1.  Images of nuclei from peripheral blood mononucleated cells stained with 4 ′,6-diamidino-2- phenylindole (DAPI) in blue, containing foci of γ -H2AX immunostained with Cyanine 3 (Cy3) in red. Figure 2.  γ -H2AX foci relative frequency found in 500 cells for four different scenarios (from left to right): at 0.5 Gy and 0.5 h, at 2 Gy and 0.5 h, at 0.5 Gy and 2 h, at 2 Gy and 2 h. The distribution of the number of observed foci varies greatly depending on the time and dose, as shown in the Fig. 2. As the dose is increased, fewer zeros appear because cells tend to have more foci. The number of foci, on the other hand, gets smaller over time. Due to the great variability over time of the foci data, the calibration curves have commonly been constructed separately for different times after e­ xposure9,10. From a practical point of view, this approach is somewhat unre- alistic given the continuity of time and the inability to accurately determine the exact time point of irradiation. Therefore, it seems a noteworthy idea to include the time as a variable in a calibration model and to consider a three-dimensional response surface, i.e. simultaneously dose and time dependent, what was recently proposed by López et al.11. Foci γ -H2AX data are ­overdispersed3,9, which means that the variance in the observed samples is higher than its mean. Due to this fact, the Poisson model, commonly used in the analysis of count data, may not be appropriate. Instead, many distributions were already proposed, in particular, negative binomial, zero-inflated Poisson and zero-inflated negative binomial models provide good results for γ -H2AX ­foci12. γ -H2AX is analyzed in PBMCs, that is composed by a mixture of various types of leukocytes, but different subsets of blood cells could have different levels of γ -H2AX foci. It has been suggested that the subset of CD4+ leukocytes show 1.5 times higher level of phosphorylation than CD19+13. Although this statement requires more investigation because more evidence is necessary, it can be a plausible explanation for the variability in the distribution of foci, hence we opted to use a mixed Poisson model in this paper. Mixture models are a useful tool for modelling data that are highly diverse, and it is thought that this behaviour is due to underlying sub-populations14. These models assume that it is not known to which subgroup a particular observation belongs, so a mixture of a few Poisson distributions with different proportions is used to describe all the data. Scientific Reports | (2022) 12:19877 | https://doi.org/10.1038/s41598-022-24331-1 2 Vol:.(1234567890) www.nature.com/scientificreports/ Figure 3.  Scheme of the whole dose estimation procedure. Moreover, many unknowns may arise in the dose estimation process, such as the impossibility of determining the exact time of irradiation or the individual response to radiation. In this scenario, the Bayesian framework may be more appropriate because it can handle many levels of u ­ ncertainty15. The Bayesian approach also allows current information to be incorporated into the inferential process, as may be the case in biodosimetry, for exam- ple when laboratories prepare a calibration curve (or surface) independent of emergency data. Therefore, some information about the distribution of the model parameters may be provided in advance of a possible accident. This study aims to investigate the use of biomarker γ -H2AX in the estimation of the radiation dose absorbed by an individual taking into account the uncertainty on their exposure time. In the second section, we introduce statistical strategies for dealing with the issues outlined before. In the third section, we apply these methods to γ -H2AX real data and we present the results. In the last part, we discuss some limitations of the statistical propos- als and we discuss some directions for further research. Statistical models In this section, we describe the proposed statistical methodology for analysing biodosimetric foci-data. The pro- cess is divided into two parts: calibration and estimation. The calibration procedure is based on sample laboratory data, implying that the irradiation exposure occurs under monitored conditions. Calibration data were obtained irradiating N peripheral mononucleated cells from one donor with radiation doses ranging from 0 to 3 Gy, and γ -H2AX foci were detected microscopically using a semi-automatic method, at different post-irradiation times from 0.5 to 24 h. The number of γ -H2AX foci, yi, was recorded, obtaining a set of data y = (y1 , . . . , yN ), assuming conditional independence among observations. The data can be found in the Supplementary material of the paper by López et al.11. These data, which serve as a type of reference for how cells behave to different doses and time points, are fitted with our statistical model to create a calibration surface. Then, blood samples from a potentially newly irradiated patient are explored, and the mean and variance of the number of foci counts in n cells are recorded. In the estimation part of the process, they are used to estimate the doses supposedly received by those irradiated persons. Calibration is only required once, and the same surface response can be used to analyse many patient samples. In fact, different laboratories can perform calibration and estimation; Fig. 3 describes the process. Calibration. As explained in the introduction, we propose a Poisson mixture model with K components to describe the sampling distribution of the number of γ -H2AX foci per cell. We assume that observations y = (y1 , . . . , yN ) are conditionally independent and have been generated by the following finite mixture  K (yi | ω, i ) ∼ ωk f (yi | ki ), (1) k=1 yi e−ki  where ω = (ω1 , . . . , ωK ), i = (1i , . . . , Ki ) and f (yi | ki ) = yi ! ki indicates the conditional Poisson prob- ability of observing yi foci givenki . The parameters 0 < ωk < 1 represent the weight of each component of the mixture and they sum up to 1, Kk=1 ωk = 1. In general, the parameter  of a Poisson random variable can be modelled in many ways using different link functions. As mentioned earlier, time since radiation exposure has a substantial influence on the presence of the number of γ -H2AX foci. For this reason, and in line with the model proposed b ­ y11,  could be defined as a function of two variables, dose and time (represented by d and t respectively),  = (d, t) = c · t u + a · t v · d, Scientific Reports | (2022) 12:19877 | https://doi.org/10.1038/s41598-022-24331-1 3 Vol.:(0123456789) www.nature.com/scientificreports/ where (a, c, u, v) are parameters. Therefore,  is a three-dimensional surface, although it is important to note that if time is fixed,  is linear with respect to d, which is consistent with previous fi ­ ndings9. For the mixture model with K components in (1), the parameters ki are defined separately for each Poisson u v component k ∈ {1, . . . , K}, ki (d, t) = ck · ti k + ak · ti k · di for i = 1, . . . , N . Since the calibration phase uses laboratory data, for a given observed number of foci yi, ti (time) and di (dose) are known. From now on to simplify the notation, θ k = (ak , ck , uk , vk ) stands for the set of parameters for component k and θ = (θ1 , . . . , θK ) for all parameters in the model. In order to find the calibration surface, these parameters, together with the weights of the mixture ω = {ω1 , . . . , ωK }, should be estimated, what we propose to do within Bayesian framework. In the Bayesian schema, all previous information about the quantities of interest are used to elicit a prior den- sity distribution π(ω, θ) for the parameters (ω, θ ). Next, the data y are registered and the likelihood L(y | ω, θ) of (ω, θ ) for the data are constructed. In our case, given the conditional independence of the observations yi , . . . , yN , the likelihood function remains, N   K u v L(y | ω, θ) = ωk f (yi | ck · ti k + ak · ti k · di ). i=1 k=1 Then the Bayes theorem provides the posterior density distribution π(ω, θ | y): π(ω, θ | y) ∝L(y | ω, θ )π(ω, θ ). As stated before, Bayesian modelling requires the specification of a prior distribution π(ω, θ ). We assume prior independence between ω and θ and, consequently π(ω, θ ) = π(ω)π(θ ). We propose to consider non-informative uniform prior distributions for θ , as they will produce a minimal influence on the inference. For ω , we decided to use a non-informative symmetric Dirichlet distribution called Perks’ prior (see Ref.16 and Supplementary material). The analytical derivation of the posterior distribution is intractable due to model’s complexity; however, the Bayesian framework provides numerous computational methods for approximating posterior distributions. In particular, sampling algorithms, which rely on Markov chain Monte Carlo (MCMC) methods, can give a good approximation of the posterior distribution. One such method is the Gibbs sampler, in which samples from π(ω, θ | y) are constructed from the conditional posterior distribution of each element in (ω, θ ) given the rest of them. Posterior estimates of (ω, θ ) are the subsequent posterior means, approximated by the corresponding sample means of the MCMC posterior distribution. Nonetheless, when dealing with a large number of model parameters, MCMC techniques are often extremely slow. If the absorbed radiation dose must be estimated immediately the technique of providing the results must be as quick as possible. Therefore, we propose to use the Laplace approximation method, which aims to find a Gaussian approximation to the posterior distribution, which results in a faster execution of the analysis. Assume that φ ∗ = (ω̂, θ̂) ∈ Rp is the mode of π(ω, θ | y), i.e. φ ∗ = arg max L(y | ω, θ )π(ω, θ). (ω,θ ) The Laplace approximation of the posterior π(ω, θ | y) provides a p-dimensional multivariate normal distribu- tion (details can be found in Supplementary material): π(ω, θ | y) ∼ Np (φ ∗ , � ˆ ), ω̂,θ̂ (2) ˆ is the inverse of the Hessian matrix of the model evaluated at (ω̂, θ̂ ) (the estimated variance-covariance where � ω̂,θ̂ matrix). The estimated parameters (ω̂, θ̂ ) will be called calibration coefficients. For a fixed time t the expected number of foci provided by the calibration surface remains,  K  K µ(d | t) = ωk k = ωk (ck · t uk + ak · t vk · d) = αt + βt d, (3) k=1 k=1 where,  K  K αt = ωk · ck · t uk = g1 (ω, θ), βt = ωk · ak · t vk = g2 (ω, θ). k=1 k=1 In (3), additional factors, such as age or gender, might be added as covariates if they were thought to be signifi- cant for dose estimation. The derived posterior distribution of the two time-dependent parameters (αt , βt ) is the output of the calibra- tion process. Using the multivariate delta method, it is straightforward to see that (αt , βt ) follows an approximated bivariate normal distribution, T ˆ ˙ N2 (g(ω̂, θ̂ ), ∇g · � π(αt , βt | y)∼ ω̂,θ̂ · ∇g ), (4) where ∇g is the gradient of g = (g1 , g2 ) evaluated at the estimated modal point (ω̂, θ̂ ). Estimation. In this part, we assume that there is available a set of new observations x = (x1 , . . . , xn ), com- ing from a new patient, called test data. In the event of a radiation emergency, these will be the counts of foci Scientific Reports | (2022) 12:19877 | https://doi.org/10.1038/s41598-022-24331-1 4 Vol:.(1234567890) www.nature.com/scientificreports/ detected in the exposed person’s blood cells, for which the absorbed dose must be estimated. After a fast explo- ration the laboratory records the mean number of foci in the test data x̄ and its standard deviation s. We can assume that2the expected number µ of foci of the patient, independently of time t, follows a normal distribution µ ∼ N (x̄, sn ), where n is the number of observations in the test data. This is a kind of nonparametric Bayesian estimate of µ justified in the Supplementary material. We use the dose relation described in (3) to relate the patient’s observed foci distribution with time t and dose d. Then, the dose can be directly determined as µ − αt (d | t) = βt , (5) where (d | t) means that dose is defined for a given value of t. Expression (5) establishes that, when the bivariate normal approximation (4) is used for approximating the posterior distribution for (αt , βt ), the posterior distribu- tion of d conditioned to t, π(d | t, y, x) can be approximated by the ratio of two dependent normal variables (see Ref.17 and Supplementary material). Given a prior information on the post-irradiation time, summarized with the prior distribution π(t), the joint posterior density of the time and dose of the new irradiated person can be expresses as π(d, t | y, x) = π(d | t, y, x)π(t). The marginal posterior density of the dose π(d | y, x) can be done by simulation accordingly with the fol- lowing algorithm: 1. Generate a value t ∗ of the prior distribution π(t). 2. Generate a value of (αt ∗ , βt ∗ ) from π(αt , βt | y). It can be done using Gibbs sampler or using the Laplace approximation. In the last case, (αt ∗ , βt ∗ ) is generated using 2 the bivariate normal distribution (4). 3. Generate a value µ∗ of the normal distribution N (x̄, sn ). Use expression (5) with the inputs µ∗, αt ∗ and βt ∗ for obtaining an estimated value of the dose d ∗ and record it. 4. Repeat steps 1–4 many times (at least 10,000). The simulated values d ∗ from π(d | y, x) allow to obtain the median and credibility intervals of the dose received by the patient. It is worth to mention that when the Laplace approximation is used the density π(d | y, x) can also be obtained by numerical integration of π(d, t | y, x) with respect to t, see Supplementary material for details. Materials and results Except for the three data sets chosen to test the models’ performance, we decided to use all of the data published ­by11 as calibration data. Radiation doses range from 0 to 3 Gy, with post-irradiation times ranging from 0.5 to 24 h. The Supplementary material contains all of the data used in this study as well as the R scripts used to obtain the results. Automatic microscopy was used to count the number of foci in 500 cells for each dose and time. The foci frequencies increase with increasing doses, as expected, but they decrease throughout time, reaching their lowest level after 24 h. The calibration data were fitted to mixture Poisson models with K = 2, 3, 4 and 5 components. Based on the Akaike information criterion ­(AIC18), we found that the best model was the mixture Poisson model with K = 4 components. Moreover, we chose a model where the parameter uk is equal in all four components and from now on it will be denoted as u. It means that for further analysis the mixture Poisson model with 4 components was chosen with mean given by k (d, t) = ck · t u + ak · t vk · d, k = 1, 2, 3, 4. leading to a 16 parameters model, the calibration coefficients, denoted as (ω̂, θ̂ ). Supplementary Figure S1 shows the profile of the calibration surface. Following our strategy, the calibration procedure additionally estimates the 136 values of the variance-covariance matrix � ˆ . Therefore, the laboratory that performs the calibration ω̂,θ̂ (indicated in Fig. 3 as Laboratory (1) generates a total of 152 figures. Based on this data, a different (or same) laboratory (denoted in Fig. 3 as Laboratory (2) can check a patient’s blood sample and determine the radiation dose they received. We used the test data from three data sets which weren’t included in the calibration data as an example of application to illustrate the methodology. Radiation doses were 0.75, 2, and 3 Gy, respectively, with post-irradi- ation durations of 4, 10, and 0.5 h. We assume that the laboratory that will conduct the patients’ blood tests will only provide aggregate foci data, such as the sample mean and variance, rather than raw data. We chose to check the method’s robustness by using three different prior distributions (π(t)) for each test data set. The first prior is a uniform distribution on the intervals (0.25, 0.75), (3, 5), and (8, 12) that is centred at the known post-irradiation times of 0.5, 4, and 10 h, respectively. The second and third options are non- standard symmetric beta distributions (see Supplementary material for details) with parameters α = 5, β = 5 and α = 100, β = 100, on the same intervals as the uniform priors. To estimate the posterior marginal distribution for the dose π(d | y, x) we have applied the Laplace approxi- mation and marginalising π(d, t | y, x). The analysis was performed using R and JAGS program and the codes can be found in Supplementary material. The results of the estimation are presented in Table 1 and the posterior densities π(d | y, x) can be seen in Fig. 4. The goodness of the Laplace approximation with respect to the Gibbs sampler has been checked providing the same results up to two decimal places (see Supplementary Table S1). It can be seen that the model performs well for first and third test data, where the estimated dose is very close to the real one independently of the prior distribution of the time since irradiation considered. However, Scientific Reports | (2022) 12:19877 | https://doi.org/10.1038/s41598-022-24331-1 5 Vol.:(0123456789) www.nature.com/scientificreports/ Test data Dose estimation Donor Dose Gy Time h Foci mean ± SE Prior distribution of time Time interval h Uniform Beta (5, 5) Beta (100, 100) 2 0.75 4.0 4.072 ± 0.230 (3, 5) 0.771 (0.567, 0.989) 0.773 (0.605, 0.948) 0.774 (0.622, 0.927) 1 2.00 10.0 4.036 ± 0.198 (8, 12) 1.38 (1.115, 1.662) 1.382 (1.159, 1.614) 1.383 (1.178, 1.591) 1 3.00 0.5 28.612 ± 0.525 (0.25, 0.75) 2.93 (2.044, 3.704) 2.955 (2.43, 3.434) 2.966 (2.785, 3.139) Table 1.  Test data and dose estimates (mean and credible interval) obtained by Laplace approximation method. Figure 4.  Marginal posterior densities π(d | y, x) for the test data obtained by Laplace approximation method. The prior distributions of the post-irradiation time are defined on the intervals: (3, 5), (8, 12), (0.25, 0.75) h. for second data set the dose is underestimated. Not surprisingly, the credible intervals are the narrowest for the prior beta distribution with α = 100, β = 100. Discussion The methods of biodosimetry are very useful in both small and large-scale a­ ccidents19. However, in both cases it is important to timely determine the radiation dose to those exposed, as it hopefully could lead to a prompt and effective treatment. Therefore the method of dose estimation should be quick and easy to use by people who provide help during ­emergencies20. The main advantage of γ -H2AX assay, comparing to others biomarkers, is its speed, because it does not require a long process of culturing cells and can provide results within a few hours of receiving a blood ­sample21, what makes γ -H2AX biomarker a good tool for rapid triage in case of a mass casualty e­ vent6. Our model is mostly suggested for acute exposures. For chronic exposures, we can expect that the observed foci at a given post-irradiation time will be a mixture of the remaining induced DSB from multiple moments in the accident. Furthermore, we think that the applicability of γ -H2AX in determining chronic exposures appears to be very limited due to its transient nature. In addition to the duration of the biological part of the procedure, we suggest also taking into account the speed of generating estimation results, because the more complex the model, the longer it takes to execute it. The Bayesian method proposed in this study, which is based on the Laplace approximation and numerical integration, is quick and accurate, and requires no simulations. When compared to the Gibbs sampling method, it produces nearly identical results, but the time required to estimate them is significantly longer. About 5000 iterations of the Gibbs sampling method were done, and the results took more than 39 h with these settings. Applying Laplace approximation allows for a significant reduction in the computation time needed to estimate the dose as the results of this method are immediate. The present model is based on data previously ­published11 where PBMCs were exposed at doses ranging from 0 to 3 Gy and γ -H2AX foci were evaluated up to 24 h post-irradiation. However, the model allows to include doses that are not evaluated at all post-irradiation times. Therefore, it can include low doses up to 24 h that will Scientific Reports | (2022) 12:19877 | https://doi.org/10.1038/s41598-022-24331-1 6 Vol:.(1234567890) www.nature.com/scientificreports/ not be analysed later, as well as high doses that will not be analysed in the first post-irradiation times, but later times such as 72 h. The ability to detect foci at late times after high-dose exposures for dose-assessment has been previously ­described10,22,23. In addition, the model would need to be validated with real radiation exposures, either accidental cases or by evaluating the kinetics of foci disappearance after whole body irradiation in cancer patients. It has been described that in vivo data coming from cancer patients showed some differences to that obtained in vitro. For example, foci analysis on thyroid cancer patients treated with radioiodine indicated that while there was an agreement between the in vitro and in vivo calibration from 0 to 2 h post-treatment, at later times the foci observed in vivo were lower to those observed in vitro indicating a different progression in the DNA r­ epair24. In the present study we are proposing a model to be used under specific circumstances, accidents involving low-LET radiation types (X- or γ -rays). There is evidence that high-LET IR radiation produces more complex and clustered DNA ­damage25 than low-LET IR, and while low-LET radiation-induced DSBs are resolved quickly DSBs induced by high-LET is s­ lower26,27. For this reason, the model here presented based on a Poisson mixture cannot be directly applied, and although Bayesian methods can be applied under different conditions, both the distribution of foci per cell, as well as their repair kinetics should be analysed previously to propose a dose estimation method based on the analysis of γ -H2AX for radiation exposures to high-LET. The model proposed in this article could be extended further. One may consider for example a scenario of dose gradient exposure. Assuming that gradient exposures are heterogeneous irradiations in the radiation field, it seems that a mixture model could be appropriate to handle these data. However, the number of mix- ture components K would probably greatly increase. Nevertheless, it is quite likely that there may occurs some interference between the differential radiosensitivity of the leukocyte subsets and the dose gradient effect. It is a topic that requires more investigation through additional experiments, the findings of which can undoubtedly be interesting. Data availability All data and computer code used for analysis during this study are included in Supplementary Information files. Received: 9 September 2022; Accepted: 14 November 2022 References 1. Rothkamm, K. & Horn, S. gamma-h2ax as protein biomarker for radiation exposure. Ann. Ist. Super. Sanita. 45(3), 265–71 (2009). 2. Rogakou, E. P., Pilch, D. R., Orr, A. H., Ivanova, V. S. & Bonner, W. M. DNA double-stranded breaks induce histone h2ax phos- phorylation on serine 139. J. Biol. Chem. 273(10), 5858–5868. https://​doi.​org/​10.​1074/​jbc.​273.​10.​5858 (1998). 3. Rothkamm, K. et al. Manual versus automated γ -h2ax foci analysis across five european laboratories: Can this assay be used for rapid biodosimetry in a large scale radiation accident?. Mutation Res. Genet. Toxicol. Environ. Mutagen. 756(1–2), 170–173. https://​ doi.​org/​10.​1016/j.​mrgen​tox.​2013.​04.​012 (2013). 4. Merkle, W. Statistical methods in regression and calibration analysis of chromosome aberration data. Radiat. Environ. Biophys. 21(3), 217–233. https://​doi.​org/​10.​1007/​BF013​23412 (1983). 5. Higueras, M., Puig, P., Ainsbury, E. A. & Rothkamm, K. A new inverse regression model applied to radiation biodosimetry. Proc. R. Soc. A. 471, 20140588. https://​doi.​org/​10.​1098/​rspa.​2014.​0588 (2015). 6. Moquet, J. et al. The second gamma-h2ax assay inter-comparison exercise carried out in the framework of the European biodo- simetry network (RENEB). Int. J. Radiat. Biol. 93(1), 58–64. https://​doi.​org/​10.​1080/​09553​002.​2016.​12078​22 (2017). 7. Mariotti, L. G. et al. Use of the γ -h2ax assay to investigate DNA repair dynamics following multiple radiation exposures. PLoS One 8(11), e79541. https://​doi.​org/​10.​1371/​journ​al.​pone.​00795​41 (2013). 8. Redon, C. E., Dickey, J. S., Bonner, W. M. & Sedelnikova, O. A. γ -h2ax as a biomarker of DNA damage induced by ionizing radia- tion in human peripheral blood lymphocytes and artificial skin. Adv. Space Res. 43(8), 1171–1178. https://​doi.​org/​10.​1016/j.​asr.​ 2008.​10.​011 (2009). 9. Einbeck, J. et al. A statistical framework for radiation dose estimation with uncertainty quantification from the γ -h2ax assay. PLoS One 13(11), e0207464. https://​doi.​org/​10.​1371/​journ​al.​pone.​02074​64 (2018). 10. Chaurasia, R. K. et al. Establishment and multiparametric-cytogenetic validation of 60co-gamma-ray induced, phospho-gamma- h2ax calibration curve for rapid biodosimetry and triage management during radiological emergencies. Mutation Res. Genet. Toxicol. Environ. Mutag. 866, 503354. https://​doi.​org/​10.​1016/j.​mrgen​tox.​2021.​503354 (2021). 11. López, J. S. et al. Establishment and validation of surface model for biodosimetry based on γ -h2ax foci detection. Int. J. Radiat. Biol. 98(1), 1–10. https://​doi.​org/​10.​1080/​09553​002.​2022.​19987​06 (2022). 12. Chilimoniuk, J. et al. countfitteR: Efficient selection of count distributions to assess DNA damage. Ann. Transl. Med.https://​doi.​ org/​10.​21037/​atm-​20-​6363 (2021). 13. Andrievski, A. & Wilkins, R. C. The response of gamma-h2ax in human lymphocytes and lymphocytes subsets measured in whole blood cultures. Int. J. Radiat. Biol. 85(4), 369–376. https://​doi.​org/​10.​1080/​09553​00090​27811​47 (2009). 14. McLachlan, G. J., Lee, S. X. & Rathnayake, S. I. Finite mixture models. Annu. Rev. Stat. Appl. 6, 355–378. https://​doi.​org/​10.​1146/​ annur​ev-​stati​stics-​031017-​100325 (2019). 15. Ainsbury, E. A. et al. Review of Bayesian statistical analysis methods for cytogenetic radiation biodosimetry, with a practical example. Radiat. Prot. Dosim. 162(3), 185–196. https://​doi.​org/​10.​1093/​rpd/​nct301 (2014). 16. Berger, J. O., Bernardo, J. M. & Sun, D. Overall objective priors. Bayesian Anal. 10(1), 189–221. https://​doi.​org/​10.​1214/​14-​BA915 (2015). 17. Pham-Gia, T., Turkkan, N. & Marchand, E. Density of the ratio of two normal random variables and applications. Commun. Stat. Theory Methods 35(9), 1569–1591. https://​doi.​org/​10.​1080/​03610​92060​06836​89 (2006). 18. Akaike, H. Information theory and an extension of the maximum likelihood principle. In Selected Papers of Hirotugu Akaike 199–213 (Springer, 1998). 19. International Atomic Energy Agency Cytogenetic Dosimetry. Applications in Preparedness for and Response to Radiation Emergen- cies (IAEA, Technical report, 2011). 20. Vinnikov, V. A., Ainsbury, E. A., Maznyk, N. A., Lloyd, D. C. & Rothkamm, K. Limitations associated with analysis of cytogenetic data for biological dosimetry. Radiat. Res. 174(4), 403–414. https://​doi.​org/​10.​1667/​RR2228.1 (2010). 21. Ainsbury, E. A. et al. Inter-and intra-laboratory comparison of a multibiodosimetric approach to triage in a simulated, large scale radiation emergency. Int. J. Radiat. Biol. 90(2), 193–202. https://​doi.​org/​10.​3109/​09553​002.​2014.​868616 (2014). Scientific Reports | (2022) 12:19877 | https://doi.org/10.1038/s41598-022-24331-1 7 Vol.:(0123456789) www.nature.com/scientificreports/ 22. Rothkamm, K. et al. Laboratory intercomparison on the γ -H2AX foci assay. Radiat. Res. 180(2), 149–55. https://​doi.​org/​10.​1667/​ RR3238.1 (2013). 23. Horn, S., Barnard, S. & Rothkamm, K. Gamma-H2AX-based dose estimation for whole and partial body radiation exposure. PLoS One 6(9), e25113. https://​doi.​org/​10.​1371/​journ​al.​pone.​00251​13 (2011). 24. Eberlein, U. et al. DNA damage in peripheral blood lymphocytes of thyroid cancer patients after radioiodine therapy. J. Nucl. Med. 57(2), 173–9. https://​doi.​org/​10.​2967/​jnumed.​115.​164814 (2016). 25. Mladenova, V., Mladenov, E., Stuschke, M. & Iliakis, G. DNA damage clustering after ionizing radiation and consequences in the processing of chromatin breaks. Molecules 27(5), 1540. https://​doi.​org/​10.​3390/​molec​ules2​70515​40 (2022). 26. Roobol, S. J. et al. Comparison of high- and low-LET radiation-induced DNA double-strand break processing in living cells. Int. J. Mol. Sci. 21(18), 6602. https://​doi.​org/​10.​3390/​ijms2​11866​02 (2020). 27. Oizumi, T., Ohno, R., Yamabe, S., Funayama, T. & Nakamura, A. J. Repair kinetics of DNA double strand breaks induced by simulated space radiation. Life (Basel) 10(12), 341. https://​doi.​org/​10.​3390/​life1​01203​41 (2020). Author contributions Conceptualization, D.M., P.P., C.A. and V.G.R.; methodology, D.M., P.P., C.A. and V.G.R.; Software, D.M. and V.G.R.; formal analysis, D.M.; resources, J.F.B. and M.P.C.; data curation, J.F.B. and M.P.C.; writing-original draft, D.M. and P.P.; writing-review and editing, D.M., P.P., C.A., V.G.R, J.F.B. and M.P.C. Funding This work was supported by the Consejería de Educación, Cultura y Deportes (Junta de Comunidades de Castilla- La Mancha (Spain)) [the Project MECESBAYES (SBPLY/17/180501/000491)]; Ministerio de Ciencia e Innovación (Spain) [research grants PID2019-106341GB-I00, RTI2018-096072-B-I00]; the Spanish Consejo de Seguridad Nuclear [BOE-A-2019-311]; and the Spanish State Research Agency [the Severo Ochoa and Marıa de Maeztu Program for Centers and units of Excellence in R &D (CEX2020-001084-M)]. Competing interests The authors declare no competing interests. Additional information Supplementary Information The online version contains supplementary material available at https://​doi.​org/​ 10.​1038/​s41598-​022-​24331-1. Correspondence and requests for materials should be addressed to D.M. Reprints and permissions information is available at www.nature.com/reprints. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. © The Author(s) 2022 Scientific Reports | (2022) 12:19877 | https://doi.org/10.1038/s41598-022-24331-1 8 Vol:.(1234567890)