- Research article
- Open Access
- Published:

# Bayesian model selection validates a biokinetic model for zirconium processing in humans

*BMC Systems Biology*
**volume 6**, Article number: 95 (2012)

## Abstract

### Background

In radiation protection, biokinetic models for zirconium processing are of crucial importance in dose estimation and further risk analysis for humans exposed to this radioactive substance. They provide limiting values of detrimental effects and build the basis for applications in internal dosimetry, the prediction for radioactive zirconium retention in various organs as well as retrospective dosimetry. Multi-compartmental models are the tool of choice for simulating the processing of zirconium. Although easily interpretable, determining the exact compartment structure and interaction mechanisms is generally daunting. In the context of observing the dynamics of multiple compartments, Bayesian methods provide efficient tools for model inference and selection.

### Results

We are the first to apply a Markov chain Monte Carlo approach to compute Bayes factors for the evaluation of two competing models for zirconium processing in the human body after ingestion. Based on *in vivo* measurements of human plasma and urine levels we were able to show that a recently published model is superior to the standard model of the International Commission on Radiological Protection. The Bayes factors were estimated by means of the numerically stable *thermodynamic integration* in combination with a recently developed copula-based Metropolis-Hastings sampler.

### Conclusions

In contrast to the standard model the novel model predicts lower accretion of zirconium in bones. This results in lower levels of noxious doses for exposed individuals. Moreover, the Bayesian approach allows for retrospective dose assessment, including credible intervals for the initially ingested zirconium, in a significantly more reliable fashion than previously possible. All methods presented here are readily applicable to many modeling tasks in systems biology.

## Background

Radioactive zirconium (Zr) isotopes are produced in large quantities in nuclear fission reactors; one of the most common fragments in uranium fission is the radioactive ^{95}Zr. For example, the estimated released ^{95}Zr activity of the Fukushima and Chernobyl accidents is considered to have detrimental health effects not only via irradiation, but also via the intake of edibles[1, 2]. The estimation of radiation doses is indispensable for risk analysis. This is true for occupational exposure[3] and patients undergoing diagnostic and therapeutic nuclear medicine[4] as well as for the public in general[5]. To calculate the radiation dose a mathematical model is required for quantifying the deposition of radioactivity from the incorporated radionuclide inside the human body. In internal dosimetry, this model is called biokinetic model as defined by the International Commission on Radiological Protection (ICRP) in[5]. Also, the ICRP put forward the current standard model, which we will simply denote the ICRP model. The parameters of this model were mostly derived from animal data. In order to obtain more reliable dose estimates for humans, the Helmholtz Zentrum München (HMGU) developed a new, physiologically more plausible biokinetic model[6]. It is based on the processing of non-radioactive Zr isotopes in 16 investigations with 12 healthy human subjects. The measurements were taken *in vivo* in plasma and urine up to 100 days after administration by application of the double tracer technique. Moreover, a global statistical analysis method has been developed and the uncertainty and sensitivity of the HMGU model parameters were analyzed[7, 8].

The biokinetic models mentioned above incorporate basic processes in the human physiological system[3, 5, 9, 10]. Mathematically, this is characterized as follows: All major human organs and tissues are simplified in the model structure as separate compartments that represent kinetically homogeneous amounts of radionuclides; the connections between organs and tissues are described via transfer rates, i.e. model parameters that represent the exchange rates between the individual mutually exclusive compartments. These multi-compartmental systems along with their transfer parameters describing the kinetic behavior of radionuclides in the human body are called *compartmental models*[5, 11]. Throughout this paper we use the terms biokinetic model and compartmental model interchangeably. The transfer of substances into and out of compartments is governed by the law of mass balance. Transfer parameters are frequently evaluated on the basis of experimental data obtained from laboratory animals and, to a lesser extent, human beings[10]. Although animal data is not directly comparable to human data, the former can nevertheless be very helpful as prior information.

In this publication, we address the problem of model inference and model selection. A Bayesian approach enables us to cover model and measurement uncertainties for a compartmental model based on human data, while simultaneously taking into account the prior information. The Bayesian framework provides a fully probabilistic approach[12]. It is grounded on the probability distribution of a problem specific parameter space conditioned on the given data. This specifies a measure of belief for all possible parameter values. Similarly – albeit not identical – to confidence intervals, Bayesian analyses provide credible sets for the parameters at stake, holding regions and limits of high parameter probability[13]. However, as they are intrinsically driven by prior informations, some care has to be taken in their interpretation[14].

For an overall assessment of the two competing biokinetic models for Zr, the previous model parameter uncertainty analysis[7, 8] is not sufficient, because uncertainties arising from the model structure were not taken into account. This shortcoming was addressed by our Bayesian approach. Considering the models themselves as a random variable allows to compute the probability for each of the models conditioned on given data. The ratio of the marginal likelihoods of two models, i.e. the ratio of the probability for the data to be produced by the specific model, is known as the *Bayes factor*, a quantity introduced by Jeffreys[15]. Performing model selection using Bayes factors is particularly useful when dealing with a few models only. While classical model selection approaches using statistics such as the AIC or likelihood ratio tests are based on single best parameter estimates[16], the Bayes factor takes into account all possible parameters values and thus deals with model uncertainty and avoids overfitting issues[17, 18]. Moreover, in contrast to classical tests, the Bayes factor provides evidence for either one of the (possibly non-nested) models by definition. With the introduction of Markov chain Monte Carlo (MCMC) methods[19–21] as tools for sampling from probability distributions, a remarkable increase in Bayesian analyses was noticed. However, due to very complex probability surfaces these approaches often struggle with sampling efficiency[22]. In order to avoid resulting convergence issues of the MCMC approach, we combined a technique called *thermodynamic integration* with a novel copula-based Metropolis-Hastings sampler[23]. This provides numerically stable results for the inference process.

The HMGU and ICRP models were compared based on *in vivo* plasma and urine data of 16 investigations of 12 human subjects[6] using Bayes factors. More precisely, the models were evaluated for each investigation individually and for the concatenated data of all investigations. The latter allows to infer transfer rates (including credible intervals) for an average individual. We furthermore provide an analysis based on the (i) plasma data and (ii) urine data individually. Throughout the analysis, the HMGU model turned out to be superior compared to the ICRP model with respect to covering the specific data. This means the HMGU model more precisely explains the given observations and therefore the processing of zirconium in the human body. We then used the HMGU model to analyze the accretion of zirconium in bones: not only did we observe a delayed aggregation but also to lesser extents than predicted by the ICRP model. Lastly, the Bayesian framework yielded credible bounds for retrospective dose assessment of an average individual, this is, based on the concatenated data of all 16 investigations. We provide a user-friendly estimation table for the prediction of initially ingested zirconium concentrations for *ex post* measurements. This impacts the estimation of intake and radiation dose, especially the bone dose, when aiming for personalized occupational monitoring data.

## Methods

### Biokinetic models for zirconium processing

We infer the biokinetics of zirconium as it is processed in the human body. The currently used compartmental model was recommended by the ICRP in[5, 10, 24] (Figure1A). It consists of eleven compartments and 15 transfer rates. Zirconium enters the body via the stomach compartment *y*_{9}and is processed until it reaches any of the two final compartments urine, *y*_{7}, or feces, *y*_{8}. Some of the transfer rates and compartments of the ICRP model are however physiologically questionable: The direct mass transport from the two bone compartments to the urinary bladder contents and upper large intestine compartments or the distinction between trabecular bone surface and cortical bone surface as such. In order to address these shortcomings, we have recently proposed the alternative HMGU model[6] combining the two bone compartments into one single compartment and replacing these mass flows by physiologically more plausible transfer rates (Figure1B). Altogether both models share eight transfer rates, which we denoted by *x*_{1},…,*x*_{8}. Transfers present in just one of the models have a unique index to facilitate distinction.

The dynamics of both models are described by a system of coupled linear first-order ordinary differential equations (ODEs): For each compartment *y*_{
j
} with time-dependent concentration *y*_{
j
}(*t*), the rate of change of *y*_{
j
}is given by

where${\mathcal{A}}_{{y}_{j}}^{+}$ is the set of indices of all transfer rates *x*_{
α
}flowing into compartment *y*_{
j
},${\mathcal{A}}_{{y}_{j}}^{-}$ the set of indices of all transfer rates flowing out of compartment *y*_{
j
}and${y}_{\left[{x}_{\alpha}\right]}$ is the compartment which is connected to *y*_{
j
}by the transfer rate *x*_{
α
}. For instance${\mathcal{A}}_{{y}_{5}}^{+}=\{8,10\}$,${y}_{\left[{x}_{8}\right]}={y}_{10}$ and${y}_{\left[{x}_{10}\right]}={y}_{1}$ in the HMGU model. We have *y*_{9}(0)=100*%* and therefore *y*_{j≠9}(0)=0*%*, this is, the complete amount of zirconium is initially contained in the stomach compartment. The explicit model specific ODE systems can be found in the Additional file1 sections 1.1 and 1.2.

### Experimental data

The human biokinetic data was measured in a stable tracer study executed at the Helmholtz Zentrum München (HMGU)[6, 25]. It includes 16 investigations of 12 healthy humans with ingestion of a investigation-specific amount of isotopically enriched stable zirconium. The administered amount was based on the individuals weight, aiming at a dose of 0.09mg stable tracer per kg body weight. Tracer concentrations were determined in blood plasma and urine. For the plasma data, samples were taken several times during the first day in increasing intervals, and more scarcely later on. Urine was collected completely in 12-24h periods on several days. The last samples were taken at 100d after tracer administration. Tracer concentrations were normalized to the respective tracer amount ingested in each investigation, such that the total ingested amount corresponds to 100*%* at *t*=0 in the stomach compartment *y*_{9}. For model development, the transfer compartment was taken to be identical with blood plasma, and concentrations therein were expressed as % per kg plasma. The plasma concentration was scaled by the total amount of plasma in the body to get absolute concentrations[26]. Urine data was expressed as excretion rate in % per day.

### Model likelihood

Technical limitations, such as missing *in vivo* measurement techniques for all involved compartments as well as noisy data introduce model uncertainties to biological systems[27]. Comparing different kinds of models based on single parameter estimates as done in maximum-likelihood approaches is thus very critical. For statistical model evaluation we here applied a Bayesian approach, taking into account the full parameter distribution: For each investigation *i* we assume that the data

follows the solution${\mathbf{c}}_{{\mathit{x}}^{k}}\left(t\right)$ of the differential equation given in (1) for any of the two models${\mathcal{M}}^{k}$ and some corresponding parameter vector *x*^{k}, where the model index *k*∈{*H* *I*}. Here,${\mathcal{M}}^{I}$ is the ICRP model and${\mathcal{M}}^{H}$ the HMGU model. Corresponding to the notation in Figure1A and1B,${\mathit{x}}^{I}=({x}_{1},\dots ,{x}_{8},{x}_{13},\dots ,{x}_{19})$ and${\mathit{x}}^{H}=({x}_{1},\dots ,{x}_{12})$. While${y}_{1}^{i,\xb7}$ indicates measurements in plasma, i.e. in the transfer compartment *y*_{1},${\stackrel{\u0307}{y}}_{7}^{i,\xb7}$ designates measurements of the excretion rate in the urine compartment *y*_{7}. There are${n}_{i}^{b}$ measurements in plasma and${n}_{i}^{u}$ measurements in urine for investigation *i*. Assuming furthermore that all data points contain normally distributed measurement errors, the investigation *i* and model *k* specific likelihood function has the form

where${c}_{{\mathit{x}}^{k}}^{b}\left({t}_{\alpha}\right)$ denotes the solution to Equation (1) for the transfer compartment *y*_{1}at time point *t*_{
α
}, corresponding to the measurement at$\left(\right)close="">{y}_{1}^{i,\alpha}$, for the parameter vector *x*^{k}. Furthermore,$\frac{d}{\mathrm{dt}}{c}_{{\mathit{x}}^{k}}^{u}\left({t}_{\beta}\right)$ is the derivative of the solution for the urine compartment *y*_{7}at time point *t*_{
β
}, corresponding to the measurement${\stackrel{\u0307}{y}}_{7}^{i,\beta}$, while *Φ*(·|*μ* *σ*) is the univariate probability density function of the normal distribution with mean *μ* and standard deviation *σ*.

The standard deviations for plasma,${\sigma}_{i}^{b}$, and for urine,${\sigma}_{i}^{u}$, are fitted investigation-specifically by simulated annealing[28] before starting the MCMC sampling process. They correspond to the combined strength of all deviations from the ODE, e.g. to the size of the measurement error as well as to the magnitude of the difference between the ODE solution and the data points that is due to natural internal fluctuations not considered by an ODE approach. The biological variability between the individual investigations is accounted for by the fact that we fit different${\sigma}_{i}^{b}$ and${\sigma}_{i}^{u}$ for each investigation i and thus get investigation-specific likelihoods. This leads to individual credible intervals for each parameter in each investigation in the MCMC sampling procedure later on.

The complete data (i.e. concatenated data) likelihood is given by${\mathcal{L}}_{\mathrm{ALL}}({\mathit{x}}^{k},k|\mathcal{D})=\prod _{i=1}^{16}{\mathcal{L}}_{i}({\mathit{x}}^{k},k\left|{\mathcal{D}}_{i}\right)$ for the complete data$\mathcal{D}=\{{\mathcal{D}}_{1},\dots ,{\mathcal{D}}_{16}\}$ where in all${\mathcal{L}}_{i}({\mathit{x}}^{k},k|{\mathcal{D}}_{i})$ the same fitted investigation independent${\sigma}_{i}^{b}={\sigma}^{b}$ and${\sigma}_{i}^{u}={\sigma}^{u}$ are used.

For the calculation of${\mathcal{L}}_{i}({\mathit{x}}^{k},k|{\mathcal{D}}_{i})$ Equation (1) has to be solved based on *x*^{k}. Since (1) is of first order, it can be written as

where${\mathbf{y}}_{{\mathit{x}}^{k}}\left(t\right)$ is the vector of all the compartments of model *k* and the time independent matrix *A*(*x*^{k}) represents all the interactions between these compartments, depending on the transfer rate values *x*^{k}. The analytical solution is then given by

The matrix exponential *e*^{A}(^{xk)t}was computed by eigenvalue decomposition using MATLAB’s eig function (see Additional file 1 section 1.3).

### Bayes factors

Specifying model specific, investigation independent prior distributions *p*(*x*^{k}|*k*) based on combined human/animal data yields the posterior distributions of the parameters for model *k* according to Bayes’ theorem:

where$p\left({\mathcal{D}}_{i}\right|k)$ denotes the marginal density obtained by integrating over the according parameter space *X*^{k}, i.e.

The quantity$p\left({\mathcal{D}}_{i}\right|k)$ is called data evidence and is the basis for the model selection process. For comparing two models *k* and *k*^{′}, we can view the model index as a random variable, apply Bayes’ theorem again and take ratios to cancel the denominator, yielding

The ratio of the two likelihoods in Equation (6) is known as the *Bayes factor*

We had no initial preference for any of the models and thus chose a uniform model prior. The Bayes factor in Equation (7) then coincides with the posterior odds ratio between the models *k* and *k*^{′} conditioned on the data${\mathcal{D}}_{i}$[18, 29].

The Bayes factor compares the posterior probability$p\left(k\right|{\mathcal{D}}_{i})$ that the data${\mathcal{D}}_{i}$ have arisen according to model *k* versus the probability$p\left({k}^{\prime}\right|{\mathcal{D}}_{i})=1-p(k\left|{\mathcal{D}}_{i}\right)$ that${\mathcal{D}}_{i}$ have arisen according to model *k*^{′}, in contrast to single parameter measures such as the AIC or the likelihood ratio test[16]. The models may be non-nested. Due to the evaluation of each model on its whole parameter space *X*^{k}(cf. Equation (5)), the Bayes factor naturally penalizes model complexity and thus prevents overfitting issues[30–32]. An interpretation of the likelihood ratio in the Bayes factor was given by Jeffreys[15], which can also be found in the Additional file 1 section 3. Most importantly, a Bayes factor${B}_{k,{k}^{\prime}}^{i}>3$*substantially* favors model *k*, while${B}_{k,{k}^{\prime}}^{i}>100$*decisively* favors model *k*. Clearly, for${B}_{k,{k}^{\prime}}^{i}<1$, model *k*^{′} is favored with evidence${B}_{{k}^{\prime},k}^{i}=1/{B}_{k,{k}^{\prime}}^{i}$.

### Prior information

Since the problem of radiation protection has been extensively researched, quite a large number of animal studies have been conducted. These yielded excessive prior knowledge for the Bayesian modeling approach. As the ICRP model is the recommended model used for dose estimation, there naturally exists information on the distribution types of the parameters involved in the model together with confidence intervals[7]. Four parameters have a lognormal distribution, five a triangular distribution and six have a normal distribution (see Additional file 1 section 2.3 for details). From the confidence intervals, the parameters of the normal and lognormal distributions were calculated. For the HMGU model detailed prior information is also available from previous studies[7, 8]. Here, eight parameters have a lognormal distribution and four a triangular one (see Additional file 1 section 2.3 for details). It is noteworthy that of the eight parameters shared in both models, *x*_{8} is the only one having different distributions in the ICRP and HMGU model. Due to a lack of prior knowledge of the dependency structure between the parameters, the multivariate prior distribution *p*(*x*^{k}|*k*) of model *k* was taken to be the product of the univariate prior distributions$p\left({x}_{r}^{k}\right|k)$ for each parameter${x}_{r}^{k}$, i.e.$p\left({\mathit{x}}^{k}\right|k)=\prod _{r}p({x}_{r}^{k}\left|k\right)$. Each univariate prior distribution was truncated at zero. While Bayes factors were computed *inter alia* for each investigation separately (see Results and discussion), the same prior information was applied throughout all investigations. This is reasonable since the priors contain information from a huge variety of preceding experiments.

### Thermodynamic integration

Computing the marginal of Equation (5) by numerical integration generally turns out to be a daunting task. There exist a variety of approaches to tackle this problem. Popular choices include Chib’s method, which is based on a Gibbs sampling scheme and therefore not always easily applicable[33] or the Reversible Jump MCMC algorithm proposed by Green[34], based on an across model search. Since the latter involves sampling of an additional model parameter, the sampling generally takes longer than sampling within the different models only. We therefore applied a within model sampling technique called *thermodynamic integration* (TI). It was proposed by Lartillot and Philippe[35] and Friel and Pettitt[36] and successfully applied e.g. in Xu *et al.*[37]. The method splits apart the computation into several intermediate steps by introducing an auxiliary “temperature” parameter *τ*∈[0,1] that governs the influence of the parameter likelihood. The basis of this approach is the *power posterior*, which is the usual posterior modified such that the likelihood is exponentiated by the temperature parameter and then normalized accordingly to obtain a probability density:

More precisely, the quantity of interest is the normalization constant

since it yields a way to compute the terms of the Bayes factor (cf. Equation (7)) by differentiating its logarithm

and then integrating both sides with respect to *τ*

according to Calderhead and Girolami[38]. This means that the natural logarithm of the marginal likelihood can be computed as the integral over the expectation of the logarithmized data likelihood within the model with respect to the power posterior. The parameter *τ* governs the flatness of the power posterior surface and, much as in the concept of path sampling[39], stabilizes the computation of Equation (5)[36]: choosing a discretization 0=*τ*_{1}<*τ*_{2}<…<*τ*_{N−1}<*τ*_{
N
}=1, we can compute the natural logarithm of the marginal likelihood$p\left({\mathcal{D}}_{i}\right|k)$ by numerically approximating the integral in Equation (11) by

Also, the expectation with respect to the power posterior in Equation (12) is approximated in the usual way by substituting it with the Monte Carlo estimate

where${\mathit{x}}_{\left(m\right)}^{k}$ denotes a sample drawn from${p}_{{\tau}_{n}}\left({\mathcal{D}}_{i}\right|k)$. For all our applications we chose a temperature schedule with *N*=30 steps according to${\tau}_{n}={\left(\frac{n-1}{N-1}\right)}^{5},n=1,\dots ,N$ to estimate$log\left(p\right({\mathcal{D}}_{i}\left|k\right))$ for each *k* and *i* as suggested by Calderhead and Girolami[38].

### Copula-based Monte Carlo sampling

The model, investigation, and temperature specific underlying Markov chain Monte Carlo (MCMC) samples were drawn using the recently introduced copula-based Metropolis-Hastings (MH) algorithm[23]. Copulas are constructs from probability theory for assessing and sampling from multivariate distributions. They are widely used in finance and ecology[40, 41]. For any absolutely continuous multivariate cumulative distribution function (cdf) *F*(*x*_{1},…,*x*_{
d
}) with marginal cdf’s *F*_{
i
}(*x*_{
i
}),*i*=1,…,*d*, joint density function *f*(*x*_{1},…,*x*_{
d
}) and marginal density functions *f*_{
i
}(*x*_{
i
}),*i*=1,…,*d*, we decompose

where$c\left({u}_{1},\dots ,{u}_{d}\right)$ is a unique copula density function defined on [0,1]^{d} with uniformly distributed marginals on [0,1]. This copula function covers the full dependency structure of the variables. In other words, every joint distribution function can be decomposed into the marginal behavior of its individual variables and a function covering its dependency structure[42]. The MH proposal function then generates problem specific proposals with an according dependence structure drawn from a pair copula distribution. Fitting the copula distribution was done in preruns containing 1,000,000 unthinned samples each. They were generated for each investigation and model separately. For back-and-forth conversion of the prerun samples and proposals[23], we naturally applied the according prior distributions of the models at hand. Choosing different conversion functions is possible, but affects the sampling performance. Before starting the final MCMC sampling procedure, the maximum a posteriori parameter estimates were computed by simulated annealing and used as initial MCMC sampling values. This makes a burn-in period dispensable. For thinning the Markov chains, i.e. for drawing approximately independent samples in the MCMC procedure, we applied the autocorrelation-based Effective Sample Size (ESS) proposed by Kass[43]. The ESS holds the number of samples left when the Markov chain is thinned such that two consecutive samples can be considered approximately independent. The copula-based MH approach is able to deal with the dependence structure in the high dimensional sampling space and allows for high proposal acceptance rates at simultaneously high ESS’s. Finally, all Bayes factors were computed based on 30,000 proposals of the copula-based MH algorithm at each *τ*_{
n
}throughout all applications.

## Results and discussion

In this section, we present the results of our analysis. We address the question which model is superiorly fitting the data. First, several general results, such as investigation dependency of the Bayes factor and effects of parameter correlations are shown, before turning to the results of the model selection, and their consequences for the HMGU and ICRP models.

### Investigation specificity of transfer rates

In radiation protection the transfer rates for the biokinetics of radionuclides in the human body are derived from data collected in various independent experiments[5]. We measured plasma and urine levels in 16 different investigations. This poses the question whether the models should be compared based on the complete dataset, or whether statistical evaluation should be done for each investigation individually. While the former approach results in one overall Bayes factor, the latter yields 16 investigation specific, not directly comparable Bayes factors. All investigations follow the same pulse-like time courses in the transfer compartment *y*_{1} while the excretion rates in the urine compartment *y*_{7}exhibit an exponential decay behavior (Figure2). However, zirconium tracer concentrations showed up to a 50-fold difference between maximal plasma concentrations, i.e. for investigation 10 (1.616*%*) and 6 (0.033*%*).

To test the hypothesis whether the diversity in concentration also effects transfer rates and therefore the estimated Bayes factors, we pairwise compared the posterior sample marginals of the MCMC run (corresponding to the samples of *τ*=1) for parameter *x*_{7}of the ICRP model between all investigations by means of a Kolmogorov-Smirnov test. Here *x*_{7}was chosen as it directly affects the observed plasma levels[8]. Except for one pair, all p-values were <6·10^{−8}, meaning that the chance of falsely rejecting the hypothesis of comparable marginals is negligible. Therefore, as the posterior marginal distributions are quite different, it can be deduced that the basis for the Bayes factor, the joint posterior distribution, can differ quite strongly w.r.t. the individuals. This indicated that each investigation should be treated separately. Nevertheless, in order to infer the transfer rates of an average subject, as currently used by the ICRP, the concatenated data had to be used. We therefore compared the HMGU and ICRP model based on both the concatenated data$\mathcal{D}=\{{\mathcal{D}}_{1},\dots ,{\mathcal{D}}_{16}\}$ and, in order to account for the biological diversity, the individual investigation specific datasets${\mathcal{D}}_{i}$ (*i*=1,…,16). This could also be the basis for further analysis of influence factors such as weight or gender.

### MCMC-based parameter estimation

Throughout, the analysis was based on 30,000 proposals for each of the 30 temperature levels in all 17 cases (one for each investigation and one for n). For the HMGU model the average ESS including one standard error, i.e. taken over all temperature levels and investigations, is 5832±405. In case of the ICRP model we obtained in average 5808±252 (approximately independent) samples from the Markov chains. This naturally implied high acceptance rates for both models. The sampling procedure thus captured the power posteriors very well.

From the posterior samples, we could derive new credible intervals for the parameters at hand as well as a new MAP estimate for an average subject which can be used if single parameter values are required (see Additional file 1 section 4.1). As can be seen in Figure3, the fit of the time courses covered the data, indicating that both models are in principle able to fit the data. This justifies our ODE approach with additive noise. However, from the fits alone it is not obvious which model is superior. Note that the credible intervals in Figure3 represent only the uncertainty based on the parameters, in contrast to measurement uncertainties accounted for by the${\sigma}_{i}^{b}$s and${\sigma}_{i}^{u}$s, which are not shown. Clearly, this uncertainty in the parameters is specific to the individual investigations or the complete data, respectively.

### Parameter correlations and model identifiability

The posterior probabilities of both the HMGU and ICRP model showed strong correlation between the parameters *x*_{7} and *x*_{8} throughout all investigations. The estimated Kendall’s *τ*’s based on the preruns were${\widehat{\tau}}_{\mathrm{HMGU}}=0.8027\pm 0.01$ and${\widehat{\tau}}_{\mathrm{ICRP}}=0.3452\pm 0.02$. This can be explained as follows: At time point *t*=0 the stomach compartment *y*_{9}is the only compartment with non-zero Zr concentration. It is exclusively connected to the small intestines *y*_{10} in all models. Therefore, all Zr compounds have to pass through *y*_{10}, which further on distributes them to the observed transfer compartment *y*_{1} via *x*_{7} or to the upper large intestines *y*_{5} via *x*_{8}. Aberrations in one of the parameters *x*_{7} or *x*_{8} thus have a direct effect on the amount of zirconium predicted for *y*_{1}. This affects the according posterior distributions. The same effect is found for the complete data$\mathcal{D}$ (compare pairwise scatterplots in Additional file 1 section 4.2). Despite the parameter dependencies, the posterior distributions of the ICRP and HMGU model are identifiable for all 16 investigations, this is, the investigation specific maximum a posteriori estimates are well defined and inferable (cf. Additional file 1 section 4.3).

### Bayesian model comparison

Using the concept of thermodynamic integration we compared the HMGU and the ICRP model based on (i) the concatenated data$\mathcal{D}=\{{\mathcal{D}}_{1},\dots ,{\mathcal{D}}_{16}\}$ and (ii) the individual investigation specific datasets${\mathcal{D}}_{i}$ (*i*=1,…,16). This results in a total of 17 Bayes factors. We found that all Bayes factors favored the HMGU model; in 14 out of the 17 cases even decisively (cf. Table1, second column, of this section and section 4 of Additional file 1).

In order to take a closer look at the contribution of the plasma and urine data to the above results, we computed additional Bayes factors based on the likelihoods${\mathcal{L}}_{i}^{b}({\mathit{x}}^{k},k|{\mathcal{D}}_{i})$ and${\mathcal{L}}_{i}^{u}({\mathit{x}}^{k},k|{\mathcal{D}}_{i})$ individually. Here, *i*=1,…,16,*ALL* and *k*∈{*I*,*H*}, where *I* represents the ICRP and *H* the HMGU model. The time courses already suggested better coverage of plasma data by the HMGU model (Figure3 above and section 4.4 of Additional file 1); for urine the situation is not that clear. This was confirmed by the Bayes factors (see Additional file 1 section 4 for sampling details): all 17 Bayes factors based on plasma data favored the HMGU model; in ten cases even decisively (Table1, third column). For the urine data, three investigations slightly favored the ICRP model (Table1, fourth column). In summary, all decisive Bayes factors are in favor of the HMGU model. While the HMGU model was never decisively rejected, the ICRP model was decisively rejected in the majority of cases. Hence, in depth analysis showed that the HMGU model is superior to the ICRP model with respect to zirconium processing in the human body. This not only holds investigation-specifically, but also based on the complete data. We additionally considered an extension of the HMGU model (see Additional file 1 section 1.2 and 4) which, however, did not show any significant improvements.

### Differences in radioactive ^{95}Zr retention in bone predicted by the HMGU and ICRP models

In internal exposure monitoring, biokinetic models were used to predict the organ retention or daily excretion of incorporated radionuclides[44]. With an interval of 120 days the radioactivity of ^{95}Zr possibly incorporated by occupational workers was routinely monitored by whole body counters. Depending on the intake route, the radiation dose of bone surfaces or colon was taken as regulatory limit for a decision if an individual is requested for person-specific monitoring[45]. In this monitoring procedure, the biokinetic model structure and parameters are used implicitly in the background. The organ retention function is the solution of the model in each compartment; the organ doses are directly related to the integral of radioactivity of ^{95}Zr in source organs over 50 years.

In order to compare the retention of ^{95}Zr as predicted by the ICRP and HMGU models, the 90% credible intervals for the time courses in the bone compartments were calculated based on the posterior samples. It is found that there is a significant difference between the two models (Figure4), where for the ICRP model we added the concentrations in the two bone compartments. The time courses were derived for stable isotopes of Zr and thus had to take the radioactive decay of ^{95}Zr with half-life of 64.032 d[46] into account. The decrease of retention in bone using the HMGU model consequently reduces the radiation dose estimate in bone in comparison to the ICRP bone dose which is currently used in monitoring.

### Retrospective dose assessment

Internal doses due to incorporated radionuclides have to be estimated with the help of biokinetic models based on indirect measurements, using for example bioassays for blood or urinary excretion. Normally, bioassay or in vivo data (e.g. radioactivity accumulated in skull or knee detected by a partial body counter) are measured after an accidental intake of radionuclides. Uncertainties of estimated doses are significant and have a large impact on the remediation and thus action costs. In contrast to conventional uncertainty analysis as performed in[7], our Bayesian approach naturally integrates the uncertainties of measured data and parameters simultaneously. This trait of the Bayesian approach is useful as it provides an estimate for the intake and its credible intervals.

For example, if the urinary excretion after accidental exposure is measured, we can compute credible intervals for the initial intake of radionuclide ^{95}Zr by exploiting the posterior distribution together with the linearity of the HMGU model. In order to be as general as possible we used the posterior samples based on the complete data. Given a posterior sample *x*^{H}, a measurement${\stackrel{\u0307}{y}}_{7}^{t}$ in [*μg*/*d*] for the urinary excretion rate of zirconium at time point *t* corresponds to a unique solution${\mathbf{c}}_{{\mathit{x}}^{H}}\left(t\right)$ of the HMGU ODE system. Due to the linearity of the ODE’s, the initial concentration${\mathbf{c}}_{{\mathit{x}}^{H}}\left(0\right)$ is by definition zero for all except the stomach compartment *y*_{9}. The latter reads${y}_{9}\left(0\right)={\stackrel{\u0307}{y}}_{7}^{t}\xb7100\%/{c}_{{\mathit{x}}^{H}}^{9}\left(t\right)$ where${c}_{{\mathit{x}}^{H}}^{9}\left(t\right)$ denotes the value of${\mathbf{c}}_{{\mathit{x}}^{H}}\left(t\right)$ in the stomach compartment at time point *t*. Now, assuming that for arbitrary posterior samples *x*^{H} the measurement${\stackrel{\u0307}{y}}_{7}^{t}$ is contained in the 90% credible interval of the solution${\mathbf{c}}_{{\mathit{x}}^{H}}\left(t\right)$ with initial condition *y*_{9}(0) as given above, lower and upper bounds for credible regions of the initial amount of zirconium taken in at *t*_{0}=0 emerge. These are based on the posterior samples. The estimated extrapolation factors for multiplication with a urine measurement (in [*μ* g/d]) after time *t* (in [*h*]) are contained in Table2 and yield the initial amount of zirconium contained in the stomach at *t*_{0}=0. For example, if an amount of${\stackrel{\u0307}{y}}_{7}^{2d}=50\mu $g/d was measured after two days, we find from Table2 that the 90% credible interval for the ingested amount lies between 0.029g and 0.059g. Since the above calculations are based on non-radioactive Zr isotopes, the results have to be corrected for dose assessment with respect to radioactive decay of the radionuclide in question, i.e. in many cases ^{95}Zr.

## Conclusions

We were the first to evaluate two competing biokinetic ODE models for zirconium processing in the human body after ingestion. These models were the current model recommended by the International Commission on Radiological Protection (ICRP) and a model developed by the Helmholtz Zentrum München (HMGU). The analysis was based on a Bayesian approach, i.e. individual Bayes factors for 16 investigations as well as a Bayes factor based on the concatenated dataset. In order to obtain reliable Monte Carlo sampling results, we combined the numerically stable thermodynamic integration with an efficient copula-based Metropolis-Hastings algorithm. In summary, the HMGU model was unequivocally superior with 14 of 17 Bayes factors being even decisive when compared to the well-established ICRP model. Also, when restricting the data on plasma and urine measurements only, we found that the HMGU model was clearly favored. The HMGU model thus best covers human data.

In contrast to the ICRP model, the HMGU model predicted a delayed accumulation of zirconium in bones which might be experimentally tested in animals in the future. Furthermore, we showed that the HMGU model can be applied for retrospective dose assessment, where the initially ingested amount of zirconium can be reconstructed including credible intervals from *ex post* urine measurements. This provides a simple hands-on tool that facilitates the decision if measures have to be taken in case of accidental exposure. In future applications the superior HMGU model together with its posterior samples can readily be used as the basis for dose estimation in internal dosimetry. The Bayesian framework for the compartmental model analysis developed in the present work is directly applicable to a personalized dose assessment and the uncertainty quantification if a person-specific monitoring is requested. More generally, the presented methodology is suitable for any ODE-based model selection task, such as the modeling of protein signaling, gene regulation, or drug processing[47], nowadays frequently put forward in systems biology[48, 49] or pharmacogenetics[50].

## References

- 1.
Eidgenössisches Nuklearsicherheitsinspektorat Informationsdienst: Radiologische Auswirkungen aus den kerntechnischen Unfällen in Fukushima vom 11.3.2011. 2011, Brugg, ENSI, Industriestrasse 19 5200 Brugg, Switzerland, http://www.ensi.ch/de/dossiers/fukushima-2/ensi-bericht-zu-fukushima-iv-radiologische-auswirkungen/

- 2.
United Nations Scientific Committee on the Effects of Atomic Radiation: Sources and Effects of Ionizing Radiation. 2008, United Nations Publications, New York

- 3.
ICRP: Limits for Intakes of Radionuclides by Workers Part 1. ICRP Publication 30. 1979, Pergamon Press, Ann.ICRP 8(4), Oxford

- 4.
ICRP: Radiation Dose to Patients from Radiopharmaceuticals. ICRP Publication 53. 1987, Pergamon Press, Ann. ICRP 18(1–4), Oxford

- 5.
ICRP: Age-dependent Doses to Members of the Public from Intake of Radionuclides (Part 1 : Ingestion dose coefficients) ICRP Publication 56. 1989, Pergamon Press, Ann. ICRP 20(2), Oxford

- 6.
Greiter M, Giussani A, Höllriegl V, Li W, Oeh U: Human biokinetic data and a new compartmental model of zirconium – A tracer study with enriched stable isotopes. Sci Total Environ. 2011, 409: 3701-3710. 10.1016/j.scitotenv.2011.06.031.

- 7.
Li W, Greiter M, Oeh U, Hoeschen C: Reliability of a new biokinetic model of zirconium in internal dosimetry Part I , Parameter uncertainty analysis. Health Phys. 2011, 101 (6): 660-676. 10.1097/HP.0b013e3181fbfba9.

- 8.
Li W, Greiter M, Oeh U, Hoeschen C: Reliability of a new biokinetic model of zirconium in internal dosimetry Part II , Parameter sensitivity analysis. Health Phys. 2011, 101 (6): 676-692.

- 9.
Guyton A, Hall J: Textbook of Medical Physiology (11th ed.). 2006, Elsevier Saunders, Philadelphia

- 10.
ICRP: Report on the Task Group on Reference Man. ICRP Publication 23. 1975, Pergamon Press, Oxford

- 11.
Jacquez J: Compartmental analysis in biology and medicine (3rd ed.). 1996, MI: BioMedware, Ann Arbor

- 12.
Clyde C, George E: Model Uncertainty. Stat Sci. 2004, 19: 81-94. 10.1214/088342304000000035.

- 13.
Marin J, Robert C: Bayesian core: a practical approach to computational Bayesian statistics. 2007, Springer Verlag, New York

- 14.
Bland J, Altman D: Bayesians and frequentists. BMJ. 1998, 317: 1151-1160.

- 15.
Jeffreys H: Some tests of significance, treated by the theory of probability. Proc Camb Philol Soc. 1935, 31: 203-222. 10.1017/S030500410001330X.

- 16.
Davison A: Statistical Models. 2003, Cambridge University Press, Cambridge

- 17.
Aris-Brosou S: How Bayes tests of molecular phylogenies compare with frequentist approaches. Bioinformatics. 2003, 19: 618-624. 10.1093/bioinformatics/btg065.

- 18.
Kass R, Raftery A: Bayes factors. J Am Stat Assoc. 1995, 90: 773-795. 10.1080/01621459.1995.10476572.

- 19.
Gelfand AE, Smith AFM: Sampling-based approaches to calculating marginal densities. J Am Stat Assoc. 1990, 85: 398-409. 10.1080/01621459.1990.10476213.

- 20.
Liu J: Monte Carlo strategies in scientific computing. 2008, Springer Verlag, New York

- 21.
Robert C, Casella G: Monte Carlo statistical methods. 2004, Springer Verlag, New York

- 22.
Ramsay J, Hooker G, Campbell D, Cao J: Parameter estimation for differential equations: a generalized smoothing approach. J R Stat Soc Series B Stat Methodol. 2007, 69: 741-796. 10.1111/j.1467-9868.2007.00610.x.

- 23.
Schmidl D, Czado C, Theis F: A vine copula based adaptive MCMC sampler for efficient inference of dynamical systems. Bayesian Anal. accepted

- 24.
ICRP: Age-dependent Doses to Members of the Public from Intake of Radionuclides (Part 2: Ingestion dose coefficients). ICRP Publication 67. 1993, Pergamon Press, Ann. ICRP 23(3–4), Oxford

- 25.
Greiter M, Höllriegl V, Oeh U: Method development for thermal ionization mass spectrometry in the frame of a biokinetic tracer study with enriched stable isotopes of zirconium. Int J Mass Spectrom. 2011, 304: 1-8. 10.1016/j.ijms.2011.02.013.

- 26.
Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P: Molecular biology of the cell (4th ed.). 2002, Garland Science, New York

- 27.
Raue A, Kreutz C, Maiwald T, Bachmann J, Schilling M, Klingmüller U, Timmer J: Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics. 2009, 25: 1923-1929. 10.1093/bioinformatics/btp358.

- 28.
Kirkpatrick S, Gelatt C, Vecchi M: Optimization by simulated annealing. Science. 1983, 220: 671-680. 10.1126/science.220.4598.671.

- 29.
Chib S, Jeliazkov I: Marginal likelihood from the Metropolis-Hastings output. J Am Stat Assoc. 2001, 96: 270-281. 10.1198/016214501750332848.

- 30.
Lodewyckx T, Kim W, Lee M, Tuerlinckx F, Kuppens P, Wagenmakers E: A tutorial on Bayes factor estimation with the product space method. J Math Psychol. 2011, 55: 331-347. 10.1016/j.jmp.2011.06.001.

- 31.
Myung I, Pitt M: Applying Occam’s razor in modeling cognition: A Bayesian approach. Psych Bull Rev. 1997, 4: 79-95. 10.3758/BF03210778.

- 32.
Pitt M, Myung I, Zhang S: Toward a method of selecting among computational models of cognition. Psychol Rev. 2002, 109: 472-491.

- 33.
Chib S: Marginal likelihood from the Gibbs output. J Am Stat Assoc. 1995, 90: 1313-1321. 10.1080/01621459.1995.10476635.

- 34.
Green P: Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika. 1995, 82: 711-732. 10.1093/biomet/82.4.711.

- 35.
Lartillot N, Philippe H: Computing Bayes factors using thermodynamic integration. Syst Biol. 2006, 55: 195-207. 10.1080/10635150500433722.

- 36.
Friel N, Pettitt N: Marginal likelihood estimation via power posteriors. J R Stat Soc Series B Stat Methodol. 2008, 70: 589-607. 10.1111/j.1467-9868.2007.00650.x.

- 37.
Xu T, Vyshemirsky V, Gormand A, von Kriegsheim A, Girolami M, Baillie G, Ketley D, Dunlop A, Milligan G, Houslay M, Kolch W: Inferring signaling pathway topologies from multiple perturbation measurements of specific biochemical species. Sci Signal. 2010, 3: ra20-10.1126/scisignal.2000517.

- 38.
Calderhead B, Girolami M: Estimating Bayes factors via thermodynamic integration and population MCMC. Comput Stat Data Anal. 2009, 53: 4028-4045. 10.1016/j.csda.2009.07.025.

- 39.
Gelman A, Meng X: Simulating normalizing constants: from importance sampling to bridge sampling to path sampling. Stat Sci. 1998, 13: 163-185.

- 40.
Min A, Czado C: Bayesian inference for multivariate copulas using pair-copula constructions. Journal of Financial Econometrics. 2010, 8: 511-546.

- 41.
Salvadori G: Extremes in nature: an approach using copulas. 2007, Springer Verlag, New York

- 42.
Kurowicka D, Joe H: Dependence Modeling: Vine Copula Handbook. 2010, World Scientific Publishing Co. Pte. Ltd, Singapore

- 43.
Neal R: Probabilistic Inference Using Markov Chain Monte Carlo Methods. 1993, Tech. rep., University of Toronto, http://www.cs.toronto.edu/∼radford/review.abstract.html

- 44.
ICRP: Individual Monitoring for Internal Exposure of Workers. ICRP Publication 78. 1997, Pergamon Press, Ann. ICRP 27(3–4), Oxford

- 45.
Bundesministerium für Umwelt Naturschutz und Reaktorsicherheit: Richtlinie für die physikalische Strahlenschutzkontrolle zur Ermittlung der Körperdosis. Teil 2: Ermittlung der Körperdosis bei innerer Strahlenexposition (Inkorporationsüberwachung) (§§40, 41 und 42 StrlSchV). 2007, Bonn

- 46.
ICRP: Nuclear Decay Data for Dosimetric Calculations. ICRP Publication 107. 2008, Pergamon Press, Ann. ICRP 38(3), Oxford

- 47.
Krumsiek J, Pölsterl S, Wittmann D, Theis F: Odefy-From discrete to continuous models. BMC Bioinformatics. 2010, 11: 233-10.1186/1471-2105-11-233.

- 48.
Becker V, Schilling M, Bachmann J, Baumann U, Raue A, Maiwald T, Timmer J, Klingmüller U: Covering a broad dynamic range: information processing at the erythropoietin receptor. Science. 2010, 328 (5984): 1404-1408. 10.1126/science.1184913.

- 49.
Raia V, Schilling M, Böhm M, Hahn B, Kowarsch A, Raue A, Sticht C, Bohl S, Saile M, Möller P, Gretz N, Timmer J, Theis F, Lehmann WD, Lichter P U K: Dynamic mathematical modeling of IL13-induced signaling in Hodgkin and primary mediastinal B-cell lymphoma allows prediction of therapeutic targets. Cancer Res. 2011, 71 (3): 693-704. 10.1158/0008-5472.CAN-10-2987.

- 50.
Zhao W, Elie V, Roussey G, Brochard K, Niaudet P, Leroy V, Loirat C, Cochat P, Cloarec S, Garaix F, Bensman A, Fakhoury M, Jacqz-Aigrain E, André J: Population pharmacokinetics and pharmacogenetics of tacrolimus in de novo pediatric kidney transplant recipients. Clin Pharmacol Ther. 2009, 86 (6): 609-618. 10.1038/clpt.2009.210.

## Acknowledgements

The authors wish to thank Nikola Müller and Andreas Raue for careful reading of the manuscript and helpful discussions. This work was supported by the Federal Ministry of Education and Research (BMBF) in its MedSys initiative (project “SysMBO”), and the European Union within the ERC grant “LatentCauses”.

## Author information

## Additional information

### Competing interests

The authors declare that they have no competing interests.

### Authors’ contributions

DS and SH implemented the algorithms, performed the computational simulations, and interpreted the results to an equal degree. DS implemented the copula-based MH algorithm. MBG provided the data. WBL and MBG supplied the models, helped analyze the data, and curated the Bayesian prior information. FJT helped to draft the manuscript and coordinated the theoretical work. DS, SH, WBL, and MBG wrote the manuscript. All authors read and approved the final manuscript.

Daniel Schmidl, Sabine Hug contributed equally to this work.

## Electronic supplementary material

## Rights and permissions

## About this article

### Cite this article

Schmidl, D., Hug, S., Li, W.B. *et al.* Bayesian model selection validates a biokinetic model for zirconium processing in humans.
*BMC Syst Biol* **6, **95 (2012) doi:10.1186/1752-0509-6-95

#### Received

#### Accepted

#### Published

#### DOI

### Keywords

- Bayesian inference
- Model selection
- MCMC sampling
- Compartmental model
- Internal dosimetry
- Systems biology