Radiation Protection Dosimetry

Radiation Protection Dosimetry Advance Access originally published online on July 25, 2008
Radiation Protection Dosimetry 2008 131(3):308-315; doi:10.1093/rpd/ncn178

This Article Abstract FREE Full Text (PDF) OA All Versions of this Article: 131/3/308    most recent ncn178v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Add to My Personal Archive Download to citation manager Google Scholar Articles by Miller, G. Articles by Guilmette, R. Search for Related Content PubMed PubMed Citation Articles by Miller, G. Articles by Guilmette, R. Social Bookmarking          What's this?

© The Author 2008. Published by Oxford University Press. All rights reserved
The online version of this article has been published under an open access model. Users are entitled to use, reproduce, disseminate, or display the open access version of this article for non-commercial purposes provided that: the original authorship is properly and fully attributed; the Journal and Oxford University Press are attributed as the original place of publication with the correct citation details given; if an article is subsequently reproduced or disseminated not in its entirety but only in part or as a derivative work this must be clearly indicated. For commercial re-use, please contact journals.permissions@oxfordjournals.org

IMPDOS (improved dosimetry and risk assessment for plutonium-induced diseases): internal dosimetry software tools developed for the Mayak worker study

G. Miller^*, L. Bertelli and R. Guilmette

Los Alamos National Laboratory, Los Alamos, NM 87545, USA

^* Corresponding author: guthrie{at}lanl.gov

Received August 27, 2007, amended May 12, 2008, accepted June 10, 2008

	ABSTRACT

TOP ABSTRACT INTRODUCTION DESCRIPTIONS OF THE METHODS Funding REFERENCES

 

A collection of software tools developed for the Mayak worker study is described. IMPDOS is software for modelling, data analysis, and activity and dose calculations using the bioassay and postmortem data from Mayak workers provided by Southern Urals Biophysics Institute. The capabilities include: parameter fitting of data for individual cases, Bayesian dose calculations using the fit results for collections of cases with extensive data as a biokinetic prior, and database storage of results for retrieval, analysis and interpretation.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION DESCRIPTIONS OF THE METHODS Funding REFERENCES

 
IMPDOS is software for modelling, data analysis, and activity and dose calculations using the bioassay and postmortem data from Mayak workers provided by Southern Urals Biophysics Institute (SUBI). This software is a part of the IMPDOS Project 2.5, Phase III, between Los Alamos National Laboratory (LANL), Los Alamos, NM, USA, and the Southern Urals Biophysics Institute (SUBI), Ozyorsk, Russia, funded by the US Department of Energy, Office of International Health Programs, EH-51. Preliminary results obtained using this software have already been reported⁽¹⁾. This software is being used to make point determinations of biokinetic parameters for a collection of Mayak cases with extensive data. The collection of biokinetic models thus determined are to be used as a discrete-empirical-Bayes⁽²⁾ biokinetic prior for the analysis of related cases using the standard LANL Bayesian methods⁽³⁾. A complete description of current results is given in a separate publication.⁽⁴⁾

IMPDOS was constructed in two programming languages. All calculations are done in FORTRAN 77, and the interface to the input files and to the results was developed using Microsoft Visual Basic 6. IMPDOS is a versatile application for fitting plutonium biokinetic models and respiratory tract parameters and for performing intake estimates, calculations of activity in compartments and calculations of tissue doses due to occupational exposures in the Mayak facilities. The fitting of the plutonium biokinetic model parameters uses the Human Respiratory Tract Model (HRTM) from the International Commission on Radiation Protection (ICRP) Publication 66⁽⁵⁾ and either the Leggett 2005 systemic model⁽⁶⁾ or the Luciani and Polig systemic model⁽⁷⁾ as model structures and starting points. The Khokhryakov model⁽⁸⁾, which contains its own lung model, may also be used.

The IMPDOS suite also uses three interlinked Microsoft access 2003 databases with their own user interfaces.

MayakData
The MS Access database MayakData.mdb contains all the data for the study. These data were obtained from S. A. Romanov and Y. V. Zaytseva at SUBI. All information regarding intake patterns, air-concentration data, bioassay data and all relevant occupational exposure data from 63 individuals are contained in this database. Of particular importance are tables of work locations and air-concentration measurements, which allow estimation of monthly (chronic) intakes for each worker.

Mayak
The MS Access database Mayak.mdb sets up, controls and stores the results of the minimum-² fitting using different modelling approaches. These point determinations of biokinetic model parameters are used to create sets of standard-format biokinetic interpolation-table files that are stored in a standard directory structure to serve as a biokinetic priors for use in Bayesian dose calculations.

MayakCalc
The MS Access database MayakCalc.mdb sets up, controls and stores the results of the Bayesian dose calculations using, as biokinetic priors, collections of biokinetics interpolation-table files created with the Mayak.mdb database. The standard LANL methods are used, numerical integration and summation⁽⁹⁾ and Markov chain Monte Carlo⁽¹⁰⁾.

This application was prepared by the Radiological Dose Assessment team, RP-2, Radiation Protection Services, LANL, Los Alamos, NM, USA. The calculational core of this software has been continuously developed and tested over many years. It has been applied in many complex situations of internal dosimetry at LANL. Activity and dose calculations for 300 times ranging from 1 to 26 000 d post-intake have been successfully intercompared against a pre-existing analytical/numerical Laplace transform method⁽¹¹⁾ and against the ICRP results⁽¹²⁾. More than 300 different biokinetic models for 34 of the most important radionuclides and intake situations, and for both ICRP-26 and ICRP-60 tissue, weighting factors have been used in this intercomparison. Other software quality-assurance exercises will be discussed in what follows.

The software is composed of various major components, which are summarised in Table 1.

View this table: [in this window] [in a new window]   Table 1. Catalogue of software components.

 
The modules SEECAL, DLSODES, ID and UF are freely available on the web (the latter from www.lanl.gov.bayesian). Persons interested in the program GN should contact its author K. Klare (kklare{at}lanl.gov). The other modules work together in complicated ways, and because the IMPDOS interface is specifically designed for the Mayak study, a significant effort would be needed to repackage this software into a form that would be self explanatory for the general user.

	DESCRIPTIONS OF THE METHODS

TOP ABSTRACT INTRODUCTION DESCRIPTIONS OF THE METHODS Funding REFERENCES

 
BK
Forward biokinetic model calculations are done by the BK program. It sets up and solves the differential equations for a general radionuclide decay chain with box (compartment) and arrow (route) transport in the body from intake to excretion. Independent biokinetics for each member of the decay chain can be assumed. If there are N_nuclides nuclides in the decay chain and for each nuclide n there are N_boxes(n) different compartments, then the general form of the equations is as follows

(1)

where N is given by

(2)

The variables y_i are the activities of the nuclide compartment for i N and the number of transformations for i > N. The numerical method allows time-dependent rate coefficients, however, except for some research studies (unusual biokinetics, effects of chelation), the rate coefficients C_i,j were independent of time. These rate coefficients are specified by the arrows in the box-and-arrow model (movements of material for the same nuclide from compartment to compartment) and the nuclear decay data (transformations of material in the same compartment from one nuclide to another). Note that nuclear transformations require that the same compartments exist for all nuclides in the decay chain below the parent nuclide. The input files for the BK program that specify the box-and-arrow model are easily readable and easily changed, as they are in the same format as the model descriptions published by the ICRP. The number of equations N is sometimes quite large, for example, for ²²⁸Th, it is about 800; however, a solution out to 26 000 d takes only 1.9 s. The BK code makes use of the Livermore Solver for ordinary differential equations with general sparse Jacobian matrix, DLSODES, by Hindmarsh and Sherman⁽¹³⁾.

The BK-calculated number of nuclear transformations in the compartments are multiplied by a user supplied specific effective energy (SEE) matrix to calculate equivalent doses to tissues. The SEE files are generated by the program SEECAL⁽¹⁴⁾. The standard compartment naming convention used by the SEECAL program must be used by BK. The basic program has been used for over a decade at LANL and has already undergone extensive QA checking⁽¹⁵⁾. Table 2 below shows comparisons between BK (LANL), the program IMBA⁽¹⁶⁾ and the DCAL⁽¹⁷⁾ program for type M and type S ²³⁹Pu (ICRP-67 plutonium systemic model), urine and fecal excretion and tissue doses.

View this table: [in this window] [in a new window]   Table 2. Comparison between urinary excretion rate values calculated by BK, IMBA and DCAL for acute inhalation of 239Pu.

 
Table 3 shows the nine cases with the poorest agreement including type S and fecal.

View this table: [in this window] [in a new window]   Table 3. Comparison between excretion rate values calculated by BK, IMBA and DCAL for acute inhalation of 239Pu.

 
Similarly, Table 4 shows the 10 cases with the poorest agreement of tissue 50 y dose.

View this table: [in this window] [in a new window]   Table 4. Comparison between dose coefficients (50 y effective dose) calculated by BK, IMBA and DCAL for acute inhalation of 239Pu.

 
These and many other comparisons show that these codes agree at a level of a few percent.

GN
Biokinetic parameters are determined by minimum-² fitting. The minimisation algorithm is the program GN, a general non-linear least squares solver written by K. Klare (unpublished). It uses a variant of the Levenberg–Marquardt (L–M) method^(18,19) to minimise the quantity

(3)

The residual r_j is the scaled deviation of the jth measurement value from the fit. The GN algorithm minimises the sum of squares of residuals, choosing damped Gauss–Newton (G–N) or L–M trust regions. It computes finite-difference derivatives, does Jacobian updates and uses an L–M method with an escape for wild updates not found in other codes. The GN algorithm is largely a G–N method with well-determined constraints on the step size (the L–M part) to keep from stepping out of the possible solution zone. The G–N part gives quadratic convergence when possible, as in most non-tricky data problems. Klare of LANL originated GN in 1986 to solve some difficult, ill-scaled, unconstrained, non-linear optimisation problems. The method is robust and efficient and it has the advantage of not requiring derivative formulas. The method was developed using an extensive set of test problems, and it was found to outperform (by getting the answer and reducing the number of evaluations) codes up to eight times its size over a wide variety of these problems.

In using GN (as well as other non-linear fitting programs), it is necessary to specify the initial values of the parameters. Convergence is not guaranteed for all starting values. Convergence was tested by always doing two fits with different parameter starting values. In the small fraction of cases where the two results disagreed, the result with the smaller value of ² was accepted.

Uncertainty calculations were done by Monte Carlo generation of multiple data sets around the initial fit values, using lognormal multiplicative and normal additive errors

(4)

where x_j is the data value, f_j the initial fit value of the jth data point, G(s) a random number with Gaussian distribution with mean 0 and standard deviation s, S_j the normalisation uncertainty and _j^(meas) the measurement uncertainty. A large number (e.g. 100) of such data sets are generated and two methods were used to calculate the resulting variations of fit parameters: (1) repeating the non-linear fit with the new data and (2) using the linearised uncertainty or error matrix; this is much faster, but only valid for sufficiently small uncertainties

(5)

where _i is the ith parameter and r_j is the jth residual, evaluated at the fit minimum. Parameter variations are calculated using the linear equation

(6)

where r_j is the variation of the jth residual caused by data variation.

Figure 1 shows a validation test case, showing a plot of true values and simulated data having normally distributed uncertainties. The true value of the first three data points is x(1) and of the second three data points is x(1) + x(2).

View larger version (14K): [in this window] [in a new window] [Download PowerPoint slide]   Figure 1. GN validation test case. The line is the true value from which data are numerically generated assuming normally distributed uncertainty. The two parameters x(1) and x(2) and their uncertainty are then calculated using the GN program.

 
Figure 2 shows the results of a GN calculation of the two parameters x(1) and x(2).

View larger version (18K): [in this window] [in a new window] [Download PowerPoint slide]   Figure 2. GN validation test case. The scattered points represent the calculated parameter values, while the star gives the true values from which the data are generated. Since this is a linear problem, the non-linear method and linear approximation agree exactly.

 
The reason the centre of the distribution shown in Figure 2 is not centred on the true value is that the data are generated about the initial fit value and not the true value, which would in practice be unknown.

FIT
The program FIT does high level control of the fitting process, interfacing with the non-linear least squares minimisation carried out by the program GN. The program FIT sets up input files for biokinetic model calculations by the program BK and allows GN to initiate BK calculations and then reads in biokinetic calculation results (as interpolation tables) to allow the calculation of ². The program FIT can vary any of the parameters of the biokinetic model calculation in order to minimise ². The Poisson/lognormal statistical model⁽²⁰⁾ of ² is used. Model parameters can be varied individually or groups of parameters may be varied at once. When parameters are grouped, a parameter in a group is varied in proportion to f^x, where x is an exponent specific to a particular parameter and f is a positive group parameter controlling all the parameters in the group. For example, if there are two parameters p₁ (exponent 1) and p₂ (exponent –1) in group A

(7)

where p₁⁽⁰⁾ and p₂⁽⁰⁾ are the initial values of parameters 1 and 2 and 0 < f_A < is the group parameter for group A. The initial values of the f parameters are always 1.

Parameters of a single intake or up to four independent intakes may be varied.

The GN program finds unconstrained minima. Constraints are built in by non-linear transformations, and the f parameters are expressed as

(8)

where X is an unconstrained variable, – < X < . To constrain the value of f to the range 1/f_max < f < f_max, the non-linear transformation

(9)

is used.

As an example, consider Mayak case number 1 from the study of 63 cases (G. Miller et al., submitted for publication). The data consist of two late-time urine measurements, the intake from air-concentration measurements⁽²¹⁾ and autopsy tissue measurements for lung, thoracic lymph nodes, skeleton, liver and other tissues (eight data points in all). Table 5 shows the intake data.

View this table: [in this window] [in a new window]   Table 5. Example of Mayak intake data (case 1).

 
The intake data are based on measured time-average air concentrations of plutonium aerosol in various rooms of the Mayak plant where the person worked, the fractions of times in each room and standard breathing assumptions⁽¹⁴⁾. It represents a variable chronic intake of plutonium from breathing ambient air. In the modelling, the intake was represented mathematically as a succession of monthly acute intakes with a fixed time profile given by the relative sizes of monthly intakes shown in Table 5 and a variable total intake amount. The monthly intakes are assumed to be represented by independent lognormal distributions with standard deviations S = 2 (Table 5 shows the mean intake). The quantity S denotes the standard deviation (of the natural logarithm) of a lognormal distribution (geometric standard deviation = exp(S)). The total intake from Table 5 is about 4.5 MBq and has S = 1.23 (variance = sum of variances of independent intakes). Table 6 shows the bioassay data.

View this table: [in this window] [in a new window]   Table 6. Example of Mayak bioassay data (case 1).

 
The quantity SD denotes the measurement uncertainty standard deviation and S denotes the lognormal standard deviation of normalisation uncertainty (geometric standard deviation = exp(S)). These uncertainty parameters are discussed in detail by G. Miller et al. (submitted for publication) and Miller⁽²⁰⁾. Because these uncertainty parameters were not given, their determination after the fact was a major task, which constitutes a large part of the work of G. Miller et al. (submitted for publication).

The fitting parameters (five in all), assuming two intake components, one type M, are as shown in Table 7.

View this table: [in this window] [in a new window]   Table 7. HRTM fitting parameters.

 
The Leggett 2005 systemic model was used. The details of the statistical analysis are given in the work of G. Miller et al. (submitted for publication).

Figure 3 shows the distribution of the intake amounts for the two components, using the method described above of Monte Carlo generation of 100 alternate realisations of the data about the initial fit value and repeats of the non-linear fit.

View larger version (27K): [in this window] [in a new window] [Download PowerPoint slide]   Figure 3. Fitted base-10 logarithms of intake amounts in kBq for Mayak case 1. A two-component intake model is assumed, with the first component having variable intake amount, slow dissolution in lung, transfer to lymph nodes and AMAD. The second, type M component has only variable intake amount.

 
The effect of the type M intake component on effective dose is shown in Figure 4. The purpose of these fits is to determine a collection of biokinetic models for later use in Bayesian dose calculations. The dose calculations for epidemiology do not use effective dose but usually lifetime absorbed dose to tissues. The reason for calculating effective dose here is to give the reader an appreciation of the magnitude of dose using the standard measure of effective dose.

View larger version (31K): [in this window] [in a new window] [Download PowerPoint slide]   Figure 4. Distributions of 50 y effective dose for fits with and without a type M component.

 
As seen, the type M component does not change the dose distribution. The unexpected result shown in Figure 4 is the large dose uncertainty even with tissue data. The data realisation with the largest dose (of 100 trials) is shown in Table 8 for a fit without the type M component.

View this table: [in this window] [in a new window]   Table 8. Largest dose case of 100 data realisations without type M component.

 
The parameter fitting illustrated in this example is used with other test cases also having extensive data to determine the variation and uncertainty of the selected biokinetic parameters across the population represented by the test cases. The simplest alternative is to ignore the parameter uncertainty. The best-fit for each test case defines a single biokinetic model. The collection of best-fit models for all test cases represents variability of biokinetics across the population. Uncertainty can be included by doing multiple fits for multiple realisations of the data as illustrated in this example. The entire collection of biokinetic models so determined are to be used as a discrete-empirical-Bayes biokinetic prior (G. Miller, submitted for publication) for Bayesian internal dose calculations for cases from the general population represented by the test cases.

	Funding

TOP ABSTRACT INTRODUCTION DESCRIPTIONS OF THE METHODS Funding REFERENCES

 
This work was part of the United States-Russian Joint Coordinating Committee for Radiation Effects Research (JCCRER) Project 2.5 and was funded under a Cooperative Agreement with the United States Department of Energy Office of International Health Programs (HS-14), Health Safety and Security Division (HSS). Funding to pay the Open Access publication charges for this article was provided by the United States Department of Energy contract for the management and operation of Los Alamos National Laboratory.

	REFERENCES

TOP ABSTRACT INTRODUCTION DESCRIPTIONS OF THE METHODS Funding REFERENCES

 

1. Bertelli L., Guilmette R. A., Miller G., Little T., Romanov S. A., Zaytseva Y. V., Vostrotin V. Modeling plutonium biokinetics for Mayak workers—a strategy for fitting sparse data. Oeh U., Roth P., Paretzke H. G., eds. Proceedings of the 9th International Conference on Health Effects of Incorporated Radionuclides: Emphasis on Radium, Thorium, Uranium and their Daughter Products, HEIR 2004, 29 November–1 December, 2004: Germany. Germany: GSF–National Research Center for Environment and Health Neuherberg.
2. Miller G. Variability and uncertainty of biokinetic model parameters—the discrete empirical bayes approximation. Radiat Prot Dosimetry. (2008, to appear, ncn180).
3. LANL. Los Alamos center for Bayesian methods in environment, safety and health. (1998) (www.lanl.gov/bayesian).
4. Miller G., Guilmette R., Bertelli L., Little T., Romanov S. A., Zaytseva Y. V. Uncertainties in internal doses calculated for Mayak workers—A study of 63 cases. In: Radiat Prot Dosimetry (2008) 131(3):316–330.[Abstract/Free Full Text]
5. International Commission on Radiological Protection. Human respiratory tract model for radiological protection. (1994) a. Oxford: Pergamon Press. ICRP Publication 66.
6. Leggett R. W., Eckerman K. F., Khokhryakov V. F., Suslova K. G., Kranhenbuhl M. P., Miller S. C. Mayak worker study: an improved biokinetic model for reconstructing doses from internally deposited plutonium. Radiat. Res (2005) 164:111–122.[CrossRef][Web of Science][Medline]
7. Luciani A., Polig E. Verification and modification of the ICRP-67 model for plutonium dose calculation. Health Phys (2000) 78(3):303–310.[Web of Science][Medline]
8. Khokhryakov V. F., Suslova K. G., Vadim V. V., Romanov S. A., Zoya S. M., Kudryavtseva T. I., Filipy R. E., Miller S. C., Krahenbuhl M. P. The development of the plutonium lung clearance model for exposure estimation of the Mayak production association, nuclear plant workers. Health Phys (2002) 82:425–431.[CrossRef][Web of Science][Medline]
9. Miller G., Inkret W. C., Martz H. F. Internal dosimetry intake estimation using Bayesian methods. Radiat. Prot. Dosim (1999) 82:5–17.[Abstract]
10. Miller G., Inkret W. C., Martz H. F., Guilmete R. Bayesian internal dosimetry calculations using Markov chain Monte Carlo. Radiat. Prot. Dosim (2002) 98:191–198.[Abstract]
11. Bertelli L., Melo D. R., Lipstein J., Cruz-Suarez R. AIDE: internal dosimetry software. Radiat. Prot. Dosim. Advance Access published on 12 March 2008.
12. International Commission on Radiation Protection. ICRP CD1: database of dose coefficients: workers and members of the public. (2002) Elsevier Health Sciences.
13. Cristy M., Eckerman K. SEECAL: Program to calculate age-dependent specific effective energies. (1993) ORNL Report ORNL/TM-12351.
14. Hindmarsh A. C., Sherman A. H. (1987) Livermore, California: Lawrence Livermore National Laboratory. LSODE.
15. Miller G., Cheng Y. S., Traub R. J., Little T. T., Guilmette R. A. Methods used to calculate doses resulting from inhalation of capstone depleted uranium aerosols. Health Phys. (2008).
16. Birchall A., Jarvis N. S., Pease M. S., Riddell A. E., Battersby W. P. The IMBA suite: integrated modules for bioassay analysis. Radiat. Prot. Dosim (1998) 79:107–110.[Abstract]
17. Eckerman K. F., et al. User's guide to the DCAL system. (2006) Oak Ridge, TN: Oak Ridge National Laboratory. ORNL/TM-2001/190.
18. Levenberg K. A method for the solution of certain problems in least squares. Quart. Appl. Math (1944) 2:164–168.
19. Marquardt D. An algorithm for least-squares estimation of nonlinear parameters. SIAM J. Appl. Math (1963) 11:431–441.[CrossRef]
20. Miller G. Statistical modeling of Poisson-lognormal data. Radiat. Prot. Dosim (2007) 124(2):155–163.[Abstract/Free Full Text]
21. Zaytseva Y. V., Tretyakov F. D., Romanov S. A., Miller G., Bertelli L., Guilmette R. A. Use of air monitoring and experimental aerosol data for intake assessment for Mayak plutonium workers. Radiat. Prot. Dosim (2007) 127:535–539.[Abstract/Free Full Text]

CiteULike    Connotea    Del.icio.us    What's this?

This article has been cited by other articles:

				G. Miller Variability and uncertainty of biokinetic model parameters: the discrete empirical Bayes approximation Radiat Prot Dosimetry, September 1, 2008; 131(3): 394 - 398. [Abstract] [Full Text] [PDF]

This Article Abstract FREE Full Text (PDF) OA All Versions of this Article: 131/3/308    most recent ncn178v1 Alert me when this article is cited Alert me if a correction is posted Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Add to My Personal Archive Download to citation manager Google Scholar Articles by Miller, G. Articles by Guilmette, R. Search for Related Content PubMed PubMed Citation Articles by Miller, G. Articles by Guilmette, R. Social Bookmarking          What's this?

Online ISSN 1742-3406 - Print ISSN 0144-8420

Copyright © 2009 Oxford University Press

Oxford Journals Oxford University Press

• Site Map
• Privacy Policy
• Frequently Asked Questions

Other Oxford University Press sites: Oxford University Press Oxford Journals China Oxford Journals Japan American National Biography Booksellers' Information Service Children's Fiction and Poetry Children's Reference Corporate & Special Sales Dictionaries Dictionary of National Biography Digital Reference English Language Teaching Higher Education Textbooks Humanities International Education Unit Law Medicine Music Online Products Oxford English Dictionary Reference Rights and Permissions Science School Books Social Sciences Very Short Introductions World's Classics