Improving Dietary Exposure Models by Imputing Biomonitoring Data through ABC Methods

1 Introduction

Many chemicals are known to accumulate in the environment and to be very slowly eliminated from the human body such as polychlorinated biphenyls, dioxins, organochlorine pesticides, chlorobenzenes, and polybrominated diphenyl ethers. They are of great concern for the public health and are listed by the Stockholm convention as Persistent Organic Pollutants [1]. To assess the risk related to these chemicals, the use of the body burden is now considered as the best way [2]. Recently, new data allow for a direct assessment of the body burden: the biomonitoring data. However, as highlighted in Lyons et al. [3], understanding the effects on public health of exposure to environmental chemicals requires establishing relationships among events along an exposure–health evaluation–risk assessment continuum. Therefore, biomonitoring data need to be linked to external exposure to be interpreted in terms of source of exposure and for risk management purposes. External exposure is usually assessed from statistical models considering current contamination levels of the different exposure media (e.g. food, air, water, etc.) and behaviors of individuals (consumption habits, activity pattern, etc.). However, these classical exposure models do not take the accumulation in the body into account. Some reverse-dosimetry or steady-state approaches aim at taking this elimination time into account through toxicokinetic models [4]. However, they do not account for the level of the past exposures and especially its evolution. While because of the long half-life of persistent chemicals, body burden measured today also reflects the past exposure of the individuals. It is well highlighted that exposure to environmental contaminants has significantly changed over the last decades due to changes in production, in use, and in successive regulations [5] but data on past exposure are rarely available. Therefore, to be able to interpret biomonitoring data when assessing the risk related to persistent chemical, a comprehensive method that takes lifelong exposure and accumulation into account is needed. It is precisely the objective of this paper to develop a method to link the biomonitoring data with the external exposure taking the accumulation and the past exposure into account.

Different statistical points are developed and combined in this paper to propose a comprehensive approach. First, we modified the Kinetic Dietary Exposure Model (KDEM) which was introduced in the context of dietary risk assessment by Bertail et al. [6]. This dynamic model takes into account the accumulation phenomenon of a given chemical in the body and a single-compartment pharmacokinetics model governing its elimination. A simple pharmacokinetics model is suitable for many environmental contaminants that accumulate mainly in lipids (e.g. MeHg, polychlorinated biphenyls, polybrominated diphenyl ethers, dioxins, etc.). The KDEM approach was used to assess exposure to dietary methylmercury (MeHg) in the French adult female population [7]. However, these results have never been validated with measured body burdens provided by biomonitoring data [8]. Moreover, the model does not take the trend of exposure into account. That is why, a modified model is developed which integrates a scaling function λ to describe the evolution of the exposure. Data on past contamination or past exposure are rarely available but it is proposed here to use biomonitoring data to estimate the parameters of this function. Given the very high complexity of the likelihood function of the new model and the computational difficulties faced when trying to optimize it, Approximate Bayesian Computation (ABC in short, Beaumont et al. [9]) is proposed to fit the numerical model from the biomonitoring data available. This approach bypasses the exact evaluation of the likelihood function to estimate the posterior distributions by means of simulations and well-chosen summary statistics. ABC has recently become popular for statistical inference in complex models such as those considered in biosciences, in population-genetics [10], and in mathematical epidemiology [11] for instance. In order to perform model selection simultaneously with the estimation of the parameter, an original two-stage hierarchical ABC-MCMC algorithm is implemented following in the footsteps of Toni et al. [12]. Specific summary statistics, especially relevant for risk assessment, are proposed based on Cramer–Von Mises statistic.

The paper is structured as follows: Section 2 outlines the key ingredients governing the modified KDEM which accounts for the possible temporal evolution in exposure. Section 3 specifies the principles of the ABC method used to fit the new dynamic model of exposure from punctual biomonitoring data. Relevant specific summary statistics for risk assessment are proposed, and the model selection for the parametric modeling of the function λ is described. Section 4 is devoted to the application of the methodology to blood concentration in dioxins measured in a population of French fishermen families.

2 Dynamic modeling of dietary exposure

The KDEM was introduced by Bertail et al. [6] in order to describe the temporal evolution of the total body burden of a chemical present in a variety of foods. Its dynamic is ruled by two components, a marked point process (MPP) which models dietary behavior and governs the accumulation of the chemical in the human body and a linear differential equation accounting for its physiological elimination.

The accumulation due to successive dietary intakes is described by the MPP (Tn,Un)n∈N, where the Tn’s are the successive times when an intake of value Un occurs (with T0=0 by convention). The Un’s are independent from the Tn’s and the inter-intake times ΔTn=Tn−Tn−1, n≥1.

The elimination process between two successive intakes is described by a single-compartment pharmacokinetic model. The total body burden of the chemical x(t) decreases with time t according to the linear differential equation:

where γ>0 denotes the elimination rate. The parameter γ is related to the biological half-life of the chemical in the body t1/2=log(2)/γ, namely the time needed for the total body burden x to decrease by half in the absence of any additional intake.

In order to integrate the evolution of the exposure into the model, a scaling function of time λ(t) is introduced in the KDEM model. This function describes the possible evolution of the dietary exposure through time and is assumed to be of a parametric form, fully determined by a vector θ. The shape of the exposure evolution is unknown. Therefore, the function λ must be represented in a flexible manner. That is why, a linear combination of functions from a B-spline basis of degree m≥1: λ(t)=∑i=0mβiBim(t) was proposed.

Finally, exposure is described by a piecewise deterministic process X(t) of initial value x0 and evolving according to eq. (1) between intake times. The total body burden values Xn=X(Tn) at intake times Tn are defined through the recurrence equation:

3 An ABC method to fit the model to biomonitoring data

We propose to implement an ABC method [9] to fit the parameter θ of the scaling function λ by comparing the distribution of the measured body burdens from data to values simulated from KDEM.

The principle of standard ABC relies on the basic rejection algorithm: a candidate value θ∗ is generated from a prior distribution Π(θ). Based on the numerical model defined by θ∗, a simulated data set X∗ is generated. It is next reduced to a vector of summary statistics s∗. One then computes a “distance”, ρ(s∗,s), between the simulated and the observed summary statistics and compares it to a tolerance level δ: if ρ(s∗,s)≤δ, θ∗ is accepted. After a chosen number of iterations, the accepted θ∗ s form the posterior distribution of θ.

4 Application to dioxin biomonitoring data in France

4.3 ABC implementation

Summary statistics of Section 3 were calculated using K=5 age classes. This number of classes guarantees a sufficient number of individuals in each class (Table 1). Observed body burdens and their distributions by age classes are depicted in Figure 1. In risk assessment, a focus is made on the most exposed individuals. Therefore, the threshold q of the distance ρ (eq. (4)) is set to 0.5. Thereby, the 50% of the most exposed individuals of each age class are used to fit the model.

Figure 1

Observed body burdens X˜ and their distributions by age classes (K=5)

Table 1

Comparison of observed and simulated body burden (mean and 95% credibility interval) distributions for each age class

Age (years)	Observed and predicted body burden (pg\kg bw)
	Size		25th p	Mean	97.5th p
18–37	124	Obs.	1,022	4,097	8,695
		Sim.	959 [854;1,152]	4,501 [4,018;5,387]	9,072 [8,463;11,615]
38–46	123	Obs.	3,253	6,358	13,881
		Sim.	1,571 [1,396;1,905]	6231 [5515;7567]	15,408 [13,662;18,685]
47–54	130	Obs.	4,254	11,709	22,327
		Sim.	1,388 [1,220;1,697]	11,397 [9,998;13,951]	25,881 [21,048;32,601]
55–61	112	Obs.	3,539	14,916	30,932
		Sim.	1,806 [2,076;2,550]	21,760 [18,945;26,729]	32,742 [26,084;45,422]
62–75	117	Obs.	9,958	20,213	43,491
		Sim.	4,686 [4,038;5,814]	22,092 [19,110;27,325]	56,144 [42,070;68,530]

The choice of tolerance δ corresponds to a trade-off between the bias and the variance of the estimate [9] and, in practice, between computability and accuracy. The tolerance is often set to be a quantile Pδ of the empirical distribution function of the distance between observed and simulated data [9, 11, 26]. In this work, Pδ=0.05 which means that the 5% of simulated X that are closest to X˜ considering all the models together are kept.

As proposed in Wegmann et al. [26], a series of 10,000 simulations have been preliminarily performed, where values for the parameter θ are each time randomly drawn from the prior distribution Π(dθ). This calibration step allows to conveniently define the tolerance δ such that Pδ=0.05. Moreover any of these simulations for which the condition ρ(X˜,X∗)≤δ is true can be used as a starting point for the chain. This removes the issue of the large number of iterations that are necessary to first satisfy the condition ρ(X˜,X∗)≤δ. These preliminary simulations are also used to adjust the proposal range by choosing an appropriate transition kernel in order to ensure a thorough exploration of the whole parameter space.

The θ which satisfy ρ(X˜,X∗)≤δ are used as a starting point for the chains of the ABC-MCMC algorithm. Thus, there is no initial burn-in period. The convergence of the chains is checked by testing the goodness-of-fit between successive intermediate posterior distributions of θ for each parameter, using a two-sample Kolmogorov–Smirnov test [27]. This test indicates that the chains have converged after 150,000 iterations. The acceptance rate is close to 38% which ensure a thorough exploration of the parameter space. The posterior distributions of each parameter are built from the accepted θ of the 150,000 iterations of two different chains. The function λ obtained from these posterior distributions for the four models are shown in black in Figure 2.

Figure 2

Function λ obtained from the prior distributions of the parameters (p5, mean, p95) in red and from posterior distributions (p5, mean, p95) in black according to the four models

4.4 Results and model validation

For the model m=2, one reaches a rate of 41% of accepted θ (against 0% for model 1, 35% for model 3, and 24% for model 4). This model is thus considered as the best regarding the selection method proposed in Section 3. To confirm this result, the Root Mean Square Error of Prediction (RMSEP) of each model is calculated. The data set is randomly divided into three subsets respecting the proportion of individuals in each age class. Then the model is fitted from two-thirds of the data set, and the RMSEP is calculated on the remaining third using the posterior distribution of θ for each model. This procedure is repeated three times, and the mean values of the RMSEP are: 13,599 for m=2, 13,535 for m=3, and 13,452 for m=4. Since, no value of θ was accepted for m=1, no RMSEP is calculated. The RMSEP decreases when the number of parameters increases. However, the values obtained remain quite comparable, and the gain of prediction provided by an increasing number of parameters does not seem significant enough. This confirms the selection of model m=2.

The scaling functions provided by each of the four tested models are shown in black in Figure 2. Since, no value of θ was accepted for m=1, there is no scaling function associated with this model. The scaling function provided by the models m=2, m=3, and m=4 are closed regarding the mean. However, the interval around the mean is much more bigger for models m=3 and m=4. The previous selection model concluded that model m=2 was the best model. The scaling function provided by this model shows that exposure has increased since the introduction of dioxins in the 1930s to a peak in the 1970s corresponding to a seven times higher exposure than the current exposure. Then, the exposure has decreased to the current exposure values.

From the obtained scaling function, the body burdens of the individuals were predicted. Thereby, the goodness-of-fit of the model is checked by comparing the distributions of the predicted body burdens and the observed body burdens. Figure 3 shows the distributions of the observed body burdens, the body burdens predicted from KDEM without scaling function and the body burdens predicted from the corrected KDEM with the scaling function considering the mean estimation of the m=2 model parameters. This comparison shows a much better fit of the corrected KDEM model to the biomonitoring data regarding the whole population.

Figure 3

Distributions of the observed body burden (solid line), the body burdens predicted by KDEM (dotted line), and the body burdens predicted by the corrected KDEM m=2 (dashed line) for the ICAR population (from the mean estimation of the model parameters)

In order to be relevant to risk assessment, the model also needs good prediction abilities. As made for the model selection, the quality of the model prediction is checked by performing a cross-validation. Body burdens for one-third of individuals are simulated with the model fitted from the remaining two-thirds. The distributions of these simulated body burdens are summarized and compared with the observed body burdens of the corresponding individuals. Results of these comparisons are given in Table 1. A distribution is considered as well predicted if the observed percentiles are included in the credibility interval of the predicted percentiles. Table 1 shows that the mean and the high percentiles are well estimated for each age class whereas for the low percentiles intervals do not contain the observed percentiles. This can be explained by the choice of specific summary statistics that only used the 50% that represent the most exposed individuals to fit the model.

5 Discussion

This paper proposes an original method based on an ABC approach to calibrate a new dynamic exposure model using biomonitoring data. This method improves chemical exposure and risk assessment as well as risk management in many ways. First, the combination of biomarker measurements with contamination levels usually exploited separately makes it possible to refine exposure assessment. The comparison of KDEM outputs with measured body burdens highlighted that the current contamination and consumption data even used in a dynamic exposure model were not sufficient to explain the actual accumulation of dioxins in the body. Indeed, environmental contamination by dioxins in France has changed with time. Individuals were thus exposed to higher intakes in the past. Unfortunately, the value of these intakes are unknown in France since no past data are available. However, the biomonitoring data from fishermen families and the modified KDEM made it possible to estimate the scaling function parameters without past data. The obtained scaling function highlights that exposure to dioxins in France peaked in the 1970s. This is fully consistent with observed evolutions in the emission of dioxin-like PCBs in the environment [28].

In France, regarding the general population, no biomonitoring data on dioxins are available. Risk assessment for the general population is thus conducted using dietary intakes provided by the Second French Total Diet Study [21]. This risk estimate does not take into account the exposure evolution and the accumulation phenomenon. This could lead to underestimate the dioxin body burdens, as observed for the Fishermen families. The KDEM corrected by the obtained scaling function can be applied to the general population intakes to predict body burdens. Therefore, assuming the exposure evolution to be the same between the fishermen population and the general population, the risk related to dioxins can be refined. The method also makes it possible to propose new risk indicators which depend on time [29]. Indeed, as individual exposure over life is described as a process X(t), the maximum exposure over a given period of time: maxt∈[0,T]X(t) or the time spent above a reference threshold u: ∫t=0TIX(t)≥udt can be calculated. A new time dimension can therefore be introduced in risk assessment.

The proposed method also improves risk management. The weakness of biomonitoring data lies in its inability to identify exposure sources and their related part to the overall exposure. Indeed, these data alone do not provide any information on the exposure route (inhalation, ingestion, etc.) or for example on the food the most involved in the body burden. Therefore, linking dietary exposure with blood concentration data constitutes therefore a useful tool for the interpretation of biomonitoring data by the risk managers. The promoted method can be used to interpret biomonitoring data from cohort studies of the US National Health and Nutrition Examination Survey [30] or the French National survey on nutrition and health [31] for instance. This can be done for dioxins but also for other substances like polychlorinated biphenyls, polybrominated diphenyl ethers, arsenic, cadmium, mercury, and pesticides.

An ABC was proposed to fit the corrected KDEM to the biomonitoring data. This application confirms that ABC methods are particularly suitable in case of complex likelihood. Indeed, the computational difficulties faced when trying to optimize the likelihood are bypassed in this method. However, ABC methods require the choice of sufficient distances and summary statistics as highlighted by Joyce and Marjoram [14] in case of high-dimensional data in genetics. This choice can be tricky, and a wrong choice can lead to fail to fit the model [15]. The distance and the summary statistics defined here were proven able to correctly predict the body burdens of the population. They were specifically defined for risk assessment and can be used for other chemicals. An original two-stage hierarchical ABC-MCMC algorithm is proposed to perform the ABC method. By considering the model as an additional parameter, the parameter estimation for each model is performed simultaneously with the model selection. However, as highlighted by Toni et al. [12] this selection procedure may penalize/favor some models because of the choices made for the prior distribution (which is also inherent in the standard Bayesian model selection) and the tolerance. For this reason, the proposed model selection was validated by a usual criterion, the RMSEP. In our application, the RMSEP confirmed the model selection made from the hierarchical algorithm.

Using the ABC-MCMC algorithm proposed by Marjoram [15], the number of proposed values to obtain a sufficient rate of accepted θ seems reasonable since after 150,000 iterations, around 57,000 θ were accepted. Otherwise when a more complicated model of exposure is used or if it is necessary to be more restrictive regarding the tolerance δ, then it would be relevant to consider alternative algorithms. Toni et al. [12] showed that the sequential Monte Carlo method, first developed by Sisson et al. [32] and derived from a sequential importance sampling algorithm from Del Moral et al. [33], may yield much better computational performance.

Regarding dioxins, a single-compartment pharmacokinetic model, determined by a single parameter, i.e. the half-life, has been proved relevant to describe the elimination process by the human body [4, 24]. This makes it possible to link exposure with biomonitoring data from a simple model with reliable statistical properties, as described by Bertail et al. [8]. Such pharmacokinetic models are suitable for many other chemicals like polychlorinated biphenyl, polybrominated diphenyl ethers, and methyl mercury, which allows for a direct use of the method proposed in this paper with these chemicals. However, for certain other chemicals, a model which more effectively describes the anatomical, physiological, physical, and chemical phenomena involved in the absorption, distribution, metabolism, and excretion of chemicals in the human body might be necessary. Physiologically based pharmacokinetic models, providing a multi-compartment description, may be a relevant choice in these situations [34].