This work describes an analysis, using a previously established chelation model, of the bioassay data collected from a worker who received delayed chelation therapy following a plutonium-238 inhalation. The details of the case have already been described in two publications. The individual was treated with Ca-DTPA via multiple intravenous injections and then nebulizations beginning several months after the intake and continuing for four years. The exact date and circumstances of the intake are unknown. However, interviews with the worker suggested that the intake occurred via inhalation of a soluble plutonium compound. The worker provided daily urine and fecal bioassay samples throughout the chelation treatment protocol, including samples collected before, during, and after the administration of Ca-DTPA. Unlike the previous two publications presenting this case, the current analysis explicitly models the combined biokinetics of the plutonium-DTPA chelate. Using the previously established chelation model, it was possible to fit the data through optimizing only the intake (day and magnitude), solubility, and absorbed fraction of nebulized Ca-DTPA. This work supports the hypothesis that the efficacy of the delayed chelation treatment observed in this case results mainly from chelation of cell-internalized plutonium by Ca-DTPA (intracellular chelation). It also demonstrates the validity of the previously established chelation model. As the bioassay data were modified to ensure data anonymization, the calculation of the true committed effective dose was not possible. However, the treatment-induced dose inhibition (in percentage) was calculated.
INTRODUCTION
Occupational internal exposure to plutonium via inhalation is a non-negligible risk, even when adequate safety controls are in place at the workplace (1). If the plutonium intake is considered significant, medical countermeasures such as chelation treatment with solution of diethylenetriaminepentaacetic acid (DTPA) in the form of calcium or zinc trisodium salt may be applied (2–4). The purpose of administering DTPA is to remove heavy metals, such as plutonium, from the body, thus reducing the internal radiation dose (2–4). However, the administration of DTPA affects plutonium's normal biokinetics inside the body, as it enhances plutonium's rate of excretion. Therefore, the standard plutonium biokinetic models cannot be applied to model the bioassay data affected by chelation (1, 5–9).
The present work aimed to interpret the bioassay data of a male worker involved in a plutonium-238 inhalation incident and treated with DTPA, as the calcium chelate (Ca-DTPA), at the French Alternative Energies and Atomic Energy Commission (CEA), France. The case has been described in-depth elsewhere by Grémy et al. (10, 11).
Previous papers describing this case (10, 11) posited that the long-term efficacy of treatment and bioassay data could be best explained by intracellular chelation of plutonium, particularly in the liver. The objective of this paper was to model the bioassay data collected after the plutonium intake, including the measurements which were affected by chelation therapy, using a recently developed (9) and validated (1) chelation model by Dumit et al. This model incorporates intracellular chelation in the skeleton and liver. Therefore, successful modeling of the data would support the hypothesis of intracellular chelation.
MATERIALS AND METHODS
Case Description
The inhalation incident, along with a full description of the chelation treatment protocol have been described elsewhere (10, 11). Briefly, the plutonium intake was detected after a routine fecal bioassay analysis. The assumption of a wound intake was excluded, as no sharp tools were used by the worker and no wound was observed. Based on all information collected, it was then assumed that the worker had an inhalation of a soluble plutonium material months before the intake was detected (10, 11). In this study the date of intake and the material solubility were treated as variables to be determined by the data.
Bioassay Data
The worker's bioassay data are summarized in Table 1, along with the normalization uncertainties used for analysis. A total of 394 data points (nData) were used. Both urine and fecal samples were collected. However, because the worker requested maximum anonymization, Grémy et al. (10, 11) multiplied the bioassay results by a factor known only to the authors (10, 11) to further ensure anonymity. This prevents calculation of the actual doses, but allows conclusions about the relative effect of the chelation treatment.
TABLE 1
Bioassay Data Used

In Fig. 1, the bioassay collection times are explained. The times are relative to the time of the first chelation treatment, which is day 0. Thus, D35 is day 35 after the first treatment.
FIG. 1
Diagram showing the time of data collections on days 35, 36, 37, etc., after the 1st chelation treatment (1st treatment on day 0, 2nd on day 21, and 3rd on day 35.) As seen in the diagram, bioassay collected “on the day of treatment” means that the 24 h collection period ended at the time of treatment. The mBq values are urinary excretion values.

In this study the bioassay data collected before, during, and after the administration of intravenously injected and nebulized Ca-DTPA were included. Several fecal bioassay data points were not included in the analysis because the results were zero, with a zero-measurement standard deviation.
The normalization uncertainties were assigned on a dataset/type of measurement basis, where the largest reasonable value of S was used for normalization uncertainties for different types of measurements (Table 1). We assigned relatively large uncertainties which are roughly similar to what we might expect to see at Los Alamos National Laboratory (LANL). Specifically, the log of the geometric standard deviation, S, for urine bioassay measurement was set to S = 0.3 and 0.5, and for fecal measurement it was set to S = 1.5 (as shown in Table 1).
Chelation Treatment
More details about the treatment regimen are available elsewhere (10, 11). Briefly, the worker started chelation therapy with Ca-DTPA at a time delayed (months) from the plutonium intake. The worker received 40 intravenous (i.v.) injections of Ca-DTPA over nearly 21 months. About a month later, the worker received 42 inhalations/nebulization of Ca-DTPA over approximately 23 months. Intravenous and nebulized therapies started respectively at 0 and 739 days after the 1st Ca-DTPA treatment. In this work, we included all the 40 intravenous injections of Ca-DTPA, and the 42 Ca-DTPA treatments administered via nebulization. Hence, a total of 82 Ca-DTPA treatments were included in the modeling.
Software
The software package IDode (12, 13) was used in this work. IDode performs internal dosimetry calculations using models that are defined by ordinary differential equations and is able to calculate parameter uncertainties using Markov chain Monte Carlo (MCMC), as seen in recent publications (1, 9, 12–17).
The software package Activity and Internal Dose Estimates (AIDE) (18) was used in this work for validating results from IDode, when possible, including validation of the dosimetric models implemented in IDode.
Biokinetic Models
All biokinetic models were constructed in IDode software (12, 13). The biokinetic models are listed in Table 2. The intake day, intake amount, one solubility parameter of the inhalation model (sp) (20), and the nebulized Ca-DTPA amounts were allowed to vary to fit the data (see Tables 2 and 3). Hence, a total of four parameters were allowed to vary. The other solubility parameters, spt and st, of the inhalation model (20) were set to zero.
TABLE 2
Biokinetic and Dosimetric Models Used

TABLE 3
Parameter Results after MCMC Calculations

The ICRP Publication 60 (22) dosimetric system was used. The dosimetry used ICRP Publication 38 (23) nuclear decay data and absorbed fraction data (i.e., alpha absorbed fractions), as published by Cristy and Eckerman (24).
Chelation Modeling Approach
This study uses a chelation modeling approach as described and used in previous publications (1, 9, 25–29), where biokinetic models describing DTPA chelates are linked together and, via second-order kinetics, combined with the plutonium systemic model (21). Together, these models allow the description of the in vivo chelation process, and when linked with the plutonium systemic model, allow the description of the transfer of material to the various organs and excretion compartments (urine and feces). The units used in this study were grams for Ca-DTPA (the non-radioactive material) and Bq for plutonium (the radioactive material) (1, 9), given that IDode is a radiation program that assumes a radioactive decay chain, where the natural unit is activity (9).
Statistical Analysis and Markov Chain Monte Carlo Method (MCMC)
The data were described mathematically by independent likelihood functions that assumed (small) normal measurement uncertainty convoluted with lognormal normalization uncertainty [combined likelihood ! to exp(–χ2/2) (1, 30)]. The MCMC method was used, which yields a chain of model parameter values, i.e., calculates a collection of alternate possible interpretations of the data in the form of the posterior distribution of parameter values (30), given the data and the assumed priors. The chain of model parameter values is used to calculate the posterior distribution of any function of the parameters, for instance, doses. The prior on each parameter was taken to be either linear-uniform over a specified range or log-scale uniform a factor of 100 greater and smaller than the starting, minimum χ2 values (30). Self consistency between the models used and the assumed likelihood functions describing the data (30) was judged by having the chain average of χ2/nData (number of data points) approximately one. The interested reader is referred to the appendix for details about the fitting procedure, convergence criteria for MCMC, etc.
RESULTS AND DISCUSSION
Material Solubility, Intake Estimation, and Modeling of Bioassay Data
The solubility of the plutonium material in the lungs was modeled using a single parameter sp, from the ICRP Publication 66 human respiratory tract model (20), which was varied to best fit the data. Given the fact the intake day was unknown, the time of intake was also varied.
The 40 intravenous injections of Ca-DTPA were included in the modeling as previously described (1, 9). As for the 42 Ca-DTPA treatments administered via nebulization, we assumed that a fixed fraction of the nebulized Ca-DTPA was instantaneously dissolved into the blood.
Gradual absorption of the nebulized Ca-DTPA, based on a priori data, was modeled but did not result in a better fit (data not shown). Because of the very rapid lung clearance, nebulized Ca-DTPA was treated as an injection of a fraction of the total inhaled amount of Ca-DTPA. Recent work confirms that only 10%-30% of inhaled Ca-DTPA reaches the blood (11).
The MCMC results are shown in Table 3, which represent an MCMC run of a total of 32,000 chain iterations carried out in parallel on 16 threads (2,000 on each thread).
Figure 2 shows scatter plots of Ca-DTPA amount (nebulized fraction) versus intake day, sp parameter versus intake amount, and committed effective dose, E(50), versus intake amount.
FIG. 2
Panel (a) shows the Ca-DTPA amount (nebulized fraction) versus the intake day. The shaded green area shows the limits of the prior on the Ca-DTPA amount and intake day. Panel (b) shows parameter sp versus the intake amount. Panels (a) and (b) show the results for the 16,000 MCMC iterations (red squares) and also shows the result when the run was continued to 32,000 iterations (blue squares). Panel (c) shows the committed effective dose, E(50), versus intake amount (black squares). Each data point represents a possible alternative interpretation of the data.

Figure 2a shows that the Ca-DTPA amount (nebulized fraction) is not limited by the prior. As one can observe by looking at the prior box shown in Fig. 2a, the Ca-DTPA amount is determined by the data.
Describing the dataset with a single solubility parameter indicated an insoluble material, in contrast to the determination of Grémy et al. (10, 11) that the material was soluble based on interviews with the worker. In this work, the best fit was obtained assuming an insoluble material, with an sp value of 7 × 10–4 (Table 3). Figure 2b shows that sp is limited by the data (not by a prior) and the values are small, meaning an insoluble material.
Figure 2c shows that the committed effective dose, E(50), is, for the most part, tightly correlated with the intake amount.
Figures 3 and 4 show the model fit to the bioassay data (urine and feces, respectively). Because of the complexity of this dataset, the χ2 breakdown by type of bioassay data, which is shown in Table 4, is easier to understand than Figs. 3b and 4b. The largest disagreement between model and data is for the 3 pre-chelation fecal samples that standout strikingly in Fig. 4a (χ2/nData = 4.09, nData = 3). Other significant disagreement (χ2/nData > 1.5) observed in this study was for the pre-chelation urine samples (χ2/nData = 2.06, nData = 3). The model performs remarkably well in describing this dataset but may have systematic differences. Going forward, one way to investigate this would be to improve the precision of bioassay measurements (discussed below in the “Study limitations” section).
FIG. 3
Model interpretation of the 238Pu data in 24 h urine. Panels (a) and (b) show different time scales on the x-axis (“Time”). Panel (a) shows the model versus the data collected before Ca-DTPA treatment (light-gray squares), the data collected after i.v. chelation with Ca-DTPA (black squares), and the data collected on the day of chelation treatment (gray triangles). The first chelation treatment occurred on day 0 and the 3rd treatment on day 35 (D35) as shown in Figs 1 and 3a (gray arrow pointing at it). Panel (b) shows the entire urine dataset. The dark gray squares are data collected after nebulized Ca-DTPA treatment. The model curve is the posterior average and the gray shade following the model curve is the standard deviation (SD) calculated by MCMC.

FIG. 4
Model interpretation of the 238Pu data in 24 h feces. These plots show the model (which had the parameter “sp” as one of the fitted parameters—see footnote “c” in Table 2, and the results shown in Table 3) versus the fecal data. Panels (a) and (b) show different time scales on the x-axis (“Time”). Panel (a) shows the model fitting the data collected before Ca-DTPA treatment (light gray squares), the data collected on the day of chelation treatment (gray triangles), and the data collected after i.v. chelation with Ca-DTPA (black squares). Panel (b) shows the entire fecal dataset. The dark gray squares are data collected after nebulized Ca-DTPA treatment. The model curve is the posterior average and the gray shade following the model curve is the standard deviation (SD) calculated by MCMC.

TABLE 4
χ2/nData by Unique Bioassay Type (Minimum χ2 Parameter Values)

Figure 3a and b shows that the chelation model is tracking the up-down pattern of chelation treatment administration in urine. However, Figure 4a and b shows only a small chelation effect in fecal excretion. To investigate whether the chelation model was accounting for fecal excretion, an exercise was conducted. The exercise consisted of fixing the value of the sp parameter to be 1/day, corresponding to a soluble material, and varying the other parameters, intake day, intake amount, and nebulized Ca-DTPA fraction, to minimize χ2. Interestingly, the result, shown in Fig. 5 is that the fecal excretion does show an up-down behavior for a soluble material, even though this is not as good of a representation of the entire dataset (χ2/nData = 2.0, nData = 394 for sp = 1 d–1, as compared to χ2/nData = 0.893, nData = 394 when sp is allowed to vary, as seen in Table 3). Thus, assumption of a soluble material increases the χ2/nData, resulting in a poorer fit of the data.
FIG. 5
Exercise results showing the model interpretation of the 238Pu data in 24 h feces when the model has the parameter “sp” set with the value of 1 d–1.

The chelation model used was developed using wound cases datasets (and with a much lesser amount of fecal data points) (9). The chelation model has been recently validated using data from an inhalation case (1). The results of the present study show the usefulness of the chelation model (1, 9) (same model structure and parameters) for another plutonium inhalation case (1).
Dose Assessment and Efficacy of Delayed Chelation Treatment
The dose assessment results, including the committed effective dose, E(50), and the committed equivalent dose to the bone surfaces, HBS, to the lungs, HLU, and to the liver, HLI, are presented in Table 5 (see also Table 3). Please note, given the fact that the bioassay data were modified to ensure data anonymization, as requested by the worker (10, 11), the calculation of the “true” committed effective dose was not possible. The fact that all bioassay results were multiplied by a single factor [which is in between 0.25 and 2, and its true value is known only by the authors of the original works (10, 11)] to ensure anonymity does not impact the validity of results for the percentage reduction of dose. Not enough plutonium was present to meaningfully contribute to depletion of Ca-DTPA, which is to say that in this parameter regime the nonlinear chelation model containing second-order kinetics is to a very good approximation linear in intake amount. Therefore, the results presented in this work are valid regardless of the true value of the factor used to multiply the bioassay data. The interested reader is referred to the appendix of Dumit et al. (9) for more discussion of this topic.
TABLE 5
Dose Assessment Results Obtained using Minimum χ 2 Values of the Parameters Shown in Table 3

The “no Ca-DTPA” calculation (Table 5) does not include Ca-DTPA treatment in the modeling. The committed effective dose E(50), and the committed equivalent doses HBS, HLU, and HLI, are significantly larger without chelation treatment, as seen in Table 5. Thus, the delayed Ca-DTPA treatment was effective in the present case of an inhalation of an insoluble plutonium material.
Studies using animal data (31, 32) have shown decorporation of plutonium from the lungs after administration of nebulized Ca-DTPA. In this study, there is no reduction of plutonium in the lungs because the second-order kinetics chelation modeling approach is not implemented in the inhalation model used (20).
Intracellular Chelation
Previous works by Grémy et al. (10, 11, 33, 34) and Dumit et al. (1, 9, 28, 29) have shown and discussed in more details regarding the ability of Ca-DTPA to chelate intracellular plutonium in the liver and skeleton. The present work supports these past studies. Figure 6a and b illustrates the removal of plutonium from both skeleton (Fig. 6a) and liver (Fig. 6b). The calculations assumed minimum χ2 parameter values as shown in Table 3.
Previous works have argued that the rapid clearance of Pu via urine results from extracellular chelation, that is chelation of Pu in blood (28, 34). However, in the model considered here, all chelation proceeds on about the same time scale, and the extracellular chelation is much smaller than intracellular chelation.
This study shows that Ca-DTPA can reduce committed whole-body and organ doses via intracellular chelation. More studies are needed to evaluate the mechanisms involved in intracellular chelation for skeleton and liver, which is out of the scope of this work. Grémy and Miccoli (35) and Dumit et al. (36) have discussed in the literature the importance of intracellular chelation, the effect of Ca-DTPA on fecal excretion, and the modeling of the biliary pathway. The interested reader is referred to these publications for more information on these topics.
Study Limitations
There were limitations with this dataset. First, the intake date and solubility of the material are unknown, with the assumption that the intake occurred via inhalation based on interviews with the worker.
Using a single solubility parameter to model the plutonium material in the lungs is undoubtedly an oversimplification, which presents potential implications on the overall findings observed in the present study. As mentioned in the Materials and Methods section and shown in the Results and Discussion section, the single solubility parameter was treated as a variable during the fitting, and, using our assumptions about data uncertainties, provides a self-consistent description of the dataset.
As mentioned in the Results and Discussion section (under “Dose Assessment and Efficacy of Delayed Chelation Treatment”), reduction of plutonium in the lungs was not observed because the second-order kinetics chelation modeling approach is not implemented in the inhalation model used (20). Although controlled animal studies have shown decorporation of plutonium from the lungs after administration of nebulized Ca-DTPA (31, 32), such implementation is out of the scope of this study. Additionally, the quality of the fit to pre-chelation fecal data is poor even with large uncertainties.
Regarding improvement of the experimental technique in future research, in addition to keeping very detailed and precise records of the timing of bioassay collections and Ca-DTPA treatment administrations, it would be better to use longer bioassay collection intervals (particularly for fecal excretion). The collection intervals should be large enough to limit biological variability to about 10%. More precise fecal measurements would permit better quantification of intracellular chelation in the liver and drive refinement of the chelation model.
SUMMARY AND CONCLUSION
The objective of this paper was to model the bioassay data collected after a plutonium intake and treatment with both intravenously injected and nebulized Ca-DTPA. This objective was successfully achieved without revision of the underlying chelation model. The present work supports the conclusions of the previous two papers (10, 11) regarding the occurrence of intracellular chelation, and demonstrates the efficacy of delayed Ca-DTPA treatment for reducing the 238Pu internal dose in tissues/organs. Finally, it demonstrated that, for modeling purposes, administration of nebulized Ca-DTPA can be treated as an i.v. injection of ≈ 40% of the administered dose.
This study involves an extremely complex and unique dataset, and the chelation model seems to successfully describe the dataset, further demonstrating its broad applicability (1, 9). It also revealed questions meriting further investigation, particularly regarding modeling the fecal excretion pathway and the interpretation of fecal bioassay data. Collaboration between researchers from different institutions is invaluable and crucial to the scientific advancement of the internal dosimetry field.
ACKNOWLEDGMENT
The authors would like to thank the anonymous reviewers for improving this manuscript with their helpful comments and recommendations.
©2023 by Radiation Research Society. All rights of reproduction in any form reserved.
REFERENCES
Appendices
APPENDIX
More Details Regarding Numerical Analysis
The analysis was carried out using the software package (12, 13) Internal Dosimetry using models based on Ordinary Differential Equations (IDode), and more information is available from the IDode help file. Because this is an exceptionally intricate numerical situation involving second-order kinetics and a very large range of numerical values, we discuss some of the details here.
First, regarding the forward-model (FM) solver times (the forward-model solutions are obtained for specified times), for the usual situation with a single intake and no chelation these are not very important, and even a relatively small number of times (e.g., 30 to 100) are sufficient to completely specify the forward-model solutions, using interpolation to fill in between the specified times. In this case with rapid changes of excretion because of chelation this is not true and a large number of solver times were required. These consisted of:
very early times (starting at –1,000 days);
very late times (out to 20,000 days);
the times of all bioassay data;
the times of all bioassay data minus 1 day (to allow calculation of 24 h excretion without interpolation);
and times relative to each chelation treatment of 0, 0.1, 0.2, 0.4, 0.8, 1.6, and 3.2 days, giving about 900 specified solver times.
Fig. A1
Forward model calculation of excretion of urine (a) and feces (b). The 24 h excretion is calculated in two ways: “Diff” means using the difference of the interpolated solutions, and “Eqs” means using the differential equations to calculate the rate of excretion and multiplying by the excretion time interval.

The FM calculation of 24 h excretion is particularly demanding in this case. Figure A1 shows both cumulative and 24 h excretion.
As shown in Fig. A1b, for fecal excretion the 24 h increment is 7 orders of magnitude smaller than the cumulative excretion. Thus, differencing the cumulative excretion to calculate the 24 h excretion requires extraordinary precision of the cumulative amount, which argues for always using the 24 h excretion calculated using the derivative of excretion from the differential equation. However, that is inaccurate in a short time after a chelation treatment, where, because of rapid time variation, nonlinearity of the cumulative excretion is important (note the difference for urine at 8 days shown in Fig. A1). As a compromise, the differential equations were used after 300 days, where the 24 h increment as a fraction of the cumulative amount starts plunging down to very small values.
TABLE A1
Information Related to Convergence of the MCMC Run. The Run Length was 2,000 Chain Iterations Calculated Independently on 16 Threads and then Combined

It is interesting to note in Fig. A1 the slight decrease of cumulative excretion at late times, which is caused by nuclear decay. As is conventionally done, nuclear decay is compensated in the calculation of 24 h excretion using either method.
Regarding the MCMC calculation of parameter uncertainty, to guarantee convergence, a very large number of forward-model calculations are needed [see discussion by Miller (30)]. The results reported here utilize 2,000 FM calculations (MCMC chain iterations) performed independently on 16 threads with the results combined. This required about 1 h using a 16-processor Acer laptop—about 0.1 sec per FM calculation.
Convergence was judged partly by the appearance of the scatter plots, looking for evidence of initialization bias (clumping) corresponding to the 16 different chain random number seeds used on the 16 different threads. Also, the comparison of results from the first half of the chain and the full chain was used to give an indication of convergence. This is shown in the Table A1 below. The priors on the parameters were either uniform (LINEAR) or log-scale uniform (LOG). The columns labeled Value and SD show either the linear average and standard deviation (LINEAR prior) or the logarithmic average and standard deviation (LOG prior).
The four varied parameters shown in Table A1 were randomly assigned to two groups of two parameters (IDode input parameter MaxParameterPerGrp = 2). This random assignment of parameters to groups was changed in a cyclical manner. At each chain iteration the parameters in only one group were only varied with the others held fixed. As discussed by Miller (30), the selection of the random walk half intervals “delta” used in the Metropolis-Rosenbluth-Taylor MCMC algorithm is very important for efficient convergence. This was done as follows.
As discussed by Miller (30), each parameter was mapped into a corresponding “theta” parameter (the cumulative prior probability) ranging from 0 to 1. The uncertainty/variability of each theta parameter in terms of a quantity called “DeltaLarge” was estimated by initial shorter MCMC runs to approximately determine the observed standard deviation (SD) of each parameter given the data and the priors and setting DeltaLarge = 3/2 SD. DeltaLarge is therefore the Delta needed to encompass essentially the full range of uncertainty of the parameter. After determining DeltaLarge for each parameter by repeating shorter MCMC runs (1,000 iterations on each thread was used) Delta is set to be DeltaLarge × DeltaFactor, where DeltaFactor is an input parameter from 0 to 1, and using repeated short runs, the size of DeltaFactor was adjusted to achieve an approximate acceptance fraction (fraction of time that the chain moves) of ½. Once DeltaFactor is set (in this case DeltaFactor was 0.5), the MCMC run was continued out to the specified number of iterations (in this case 2,000 on each thread). With a larger number of iterations, the value of acceptance fraction changes somewhat, and in this case the final acceptance fraction was 0.42.
Thus, with given DeltaLarge for each parameter, the single input parameter DeltaFactor was used to set the acceptance fraction to be about ½. A large acceptance fraction goes along with small delta and the chain exploring only a small region of parameter space and slow convergence, while a small acceptance fraction with large delta more rapidly explores a larger region of parameter space but with a larger fraction of the time where the chain does not move.
Calling for a dose calculation doubles the number of differential equations (in this case from 64 to 128) because of needing to calculate the number of nuclear transformations in each compartment in addition to the activity in each compartment. The MCMC calculation of Dose used the parameter values of the combined saved chains (the combined “tape” file) without rerunning the chain, so the more complex FM calculation was done only for the saved tape-file parameter values.