Estimation of occupational exposure to 2,3,7,8-tetrachlorodibenzo-p-dioxin using a minimal physiologic toxicokinetic model.

In this study we investigated estimation of occupational exposure to 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) based on a minimal physiologic toxicokinetic model in humans. Our purpose was to obtain a mathematical tool for dose-response studies based on human data. We first simplified an existing model of TCDD kinetics in humans and estimated its parameters (i.e., liver elimination and background input of TCDD) using repeated measures of serum dioxin taken in Vietnam veterans (Ranch Hand data and data from an unexposed reference group). We carried out computer simulation and estimation of the model parameters both under a nonlinear weighted least-squares model (naive pooled data approach) and under a nonlinear mixed-effects model. The best parameter estimates were obtained with log-transformed data under a mixed-effects model: liver elimination parameter kf = 0.022 days-1 (95% confidence interval [CI] = 0.020, 0.024), and background input rate input = 0.1251 pg/kg/day (95% CI = 0.071, 0.179). The dioxin kinetic model and its estimated parameters were then used to provide dose estimates for a cohort of workers with exposure to TCDD at chemical plants in the United States. First, the model was used to estimate the rate of occupational intake of TCDD in a subset of the cohort consisting of 253 subjects for whom one measure of serum TCDD was available. A model of change in body-mass index over time was also identified for this subsample. The occupational exposure rate was estimated by linear regression using the above values of kinetic parameters and assuming an initial condition for serum TCDD of 7 ppt, i.e., the average level found in unexposed workers. The estimate of the occupational exposure parameter was 232.7 pg/kg/day (95% CI 192, 273). This value can be applied to the full cohort to obtain for each cohort member the time course of serum dioxin concentration from which exposure indices can be derived. Sensitivity coefficients to model parameters (background input, kf, occupational exposure, and the assumed TCDD concentration at hire) allow for a convenient recalculation of the serum TCDD curve and of the derived exposure indices for different assumed values of the model parameters.


Introduction
2,3,7,8-Tetrachlorodibenzo-p-dioxin of TCDD has been studied in recent years (TCDD) dose-response analyses based on using different mathematical modeling human data require the use of a kinetic approaches, ranging from statistical model for TCDD. The kinetic behavior regression models to comprehensive descriptions of the biologic pathways of TCDD.
Statistical regression models aim at providing a black-box description of the variables that influence TCDD kinetics. These models generally assume a fixed halflife for TCDD and additional covariates account for deviations from the fixed halflife model. One of the first attempts to describe TCDD kinetics in humans, in particular fractional clearance rates, was based on a one-compartment, time invariant, model. An estimate of 7.1 years (95% confidence interval [CI] = 5.8-9.6 years) was obtained for the half-life of TCDD in a group of 36 Ranch Hand (RH) veterans (1). This estimate changed to 11.3 years (95% CI= 10.0-14.1 years) by extending the study group to 337 veterans (2). A mixed-effects modeling approach was later adopted for describing TCDD elimination rate as a linear function of individual percent body fat, change in percent body fat, and age (3). The reported unadjusted estimated half-life is 8.7 years (95% CI = 8.0-9.5 years), although there was a statistically significant increase with increasing body fat but not with age or relative changes in body fat. Regression models incorporating one-compartment first-order kinetics for TCDD have been used also by others (4)(5)(6).
Physiologic models focus on mechanistic relations between variables. Their scope is often to provide a highly detailed description of the network of biochemical and biophysical processes related to TCDD (7)(8)(9). Animal experimental data are the source of parameter values and of model validation. The identification of physiologic models for TCDD in humans is complicated by the difficulties in obtaining the necessary observations to estimate parameter values for the modeled processes. Recently, a toxicokinetic model for TCDD based on a minimal physiologic construct has been proposed (10). The model assumes a fixed fractional clearance rate for hepatic TCDD degradation with a daily TCDD intake proportional to body weight (bw). The model accounts for variations in TCDD serum levels due to variations in body mass, even in the absence of any change in the rate of dioxin intake or elimination from the body. The model is designed to describe long-term behavior of dioxin, and it does not account for fast dynamics nor for liver sequestration and binding of TCDD (11). A diagram of the Environmental Health Perspectives * Vol 106, Supplement 2 * April 1998 THOMASETH AND SALVAN model is shown in Figure 1. The model focuses on the dynamics of TCDD in the lipid fractions, and it assumes a continuous equilibrium of TCDD concentration between liver lipids and the remaining lipid compartments. The model requires the estimation of lipid compartment volumes of adipose tissue, liver, and other tissues by means of anthropometric formulas involving body weight and height (12) and of assigned constants of fractional tissue lipid content (13). The model parameters can be estimated from repeated measures of serum TCDD.
We begin by reviewing the model presented by Dankovic et al. (10). Next, we reformulate the model as a time-variant compartmental model under the simplifying assumption of constant body height. In the reformulated model, the fractional clearance and serum concentration of TCDD depend on the individual time course of body-mass index (BMI). We then focus on estimating the model parameters from sparse data with the goal of determining the characteristic population kinetic parameters of TCDD, including their interindividual variability. We then extend the model to include occupational exposure, and we estimate the resulting additional parameter. Finally, we use the model to obtain serum TCDD profiles over time and exposure indices (area under the curve) for members of an occupational cohort (14).

Materials and Methods
Subjc Ranch Hand Data Set. Operation Ranch Hand was the unit responsible for aerial spraying of herbicides in Vietnam from 1962 to 1971 (7). Veterans are involved in a 20-year prospective study. Descriptions of the Ranch Hand data are presented by Mickalek et al. (3) and Wolfe et al. (15). The data set made available to us by Michalek contained 2362 observations on male subjects, 1008 of them pertaining to exposed (U.S. Air Force veterans of Operation Ranch Hand) and 1354 to unexposed (reference) subjects (U.S. Air Force personnel who served in Southeast Asia but who were not involved in herbicide spraying operations). From both groups we selected only observations with repeated measures of serum TCDD. We adopted the following data restrictions: For the exposed subjects' data set (RH), we used only the observations from the follow-up study, which is restricted to subjects with a 1987 TCDD serum level greater than 10 ppt (2). There are 279 observations with measurements taken in 1982, 1987, and 1992. One observation had one nonquantifiable measurement, which we discarded. One observation had a nondetectable value as the third data point, which we discarded as an outlier. This leaves 277 exposed subjects with three data points. Additionally, we included 65 RH observations with two data points, taken in 1982 and 1987.
In the reference group there are 43 observations with two data points taken in 1987 and 1992. Reference subjects with two serial TCDD measurements were not selected based on any known criteria. In particular, the availability of the second data point was not dependent on the value of the first, nor was the value of the first measurement known to individuals who volunteered for the second measurement (16). Of the 43 observations 9 were subjects with nondetectable values. After an initial attempt at imputing these 9 values from knowledge of the detection limits, we discarded the observations as their detection limit was very high in some cases, possibly implying a low precision of the measurement. One subject displayed a large variation between the two dioxin measurements and was discarded as an outlier. This left 33 observations from the reference group.
A random sample of the RH data with the 1987 serum TCDD less than 10 ppt was offered an additional TCDD measure in 1992 (16). There were 16 of these observations with detectable values. One additional observation also had the 1982 level measured. These 17 subjects provide information on background exposure input, and thus they function as unexposed subjects.
Summary statistics for the restricted data set are reported in Table 1. NIOSH Subcohort Data. Estimation of TCDD occupational exposure was carried out on a subsample of 253 male workers from the National Institute for Occupational Safety and Health (NIOSH) Fkf Figure 1. One-compartment representation of the TCDD kinetic model. The model is based on the following assumptions: a) dynamic equilibrium of TCDD concentrations between different body lipid distribution volumes; b) first order elimination proportional to TCDD liver content; and c) daily intake proportional to body weight. cohort (14,17) for whom a single measure of serum TCDD was available, usually taken long after termination of employment. Two data points for height and weight are available, at hire and at the time of the exam. The NIOSH cohort subsample of 253 workers contains 42 missing values for the measures of height and weight at hire. BMI values at hire for these subjects were imputed by the conditional means (Buck's) method (18). The method estimates the missing values by a linear regression on available predictors. The predictors we used were BMI at the time of the TCDD measure and the time interval between hire and TCDD measurement. We checked the results against the estimate of occupational exposure obtained from the complete data of 211 observations. Summary statistics for the NIOSH subcohort are reported in Table 2. NIOSH Cohort Data. Our purpose in obtaining estimates of occupational exposure to TCDD based on the minimal physiologic toxicokinetic (MPTK) model was to eventually obtain predicted serum profiles of TCDD over time and other derived exposure indices for the NIOSH cohort (14). The cohort consists of 5172 male workers employed at 12 chemical plants in the United States. For the purpose of our analysis, the following observations should be exduded: set 1, n= 202 individuals without detailed work history; set 2, n= 983 with missing height or weight. Sets 1 and 2 are not mutually exclusive, leaving a total of 4053 subjects. In this study we show only an example of application.

MPTK Modeling ofTCDD in Humans
The MPTK model of TCDD in humans proposed by Dankovic et al. (10) provides a concise description of time variations in serum lipids concentration of TCDD in terms of liver degradation and variations in body lipids. The scope of the model is to describe long-term kinetics of TCDD in lipid fractions of several tissues by a minimal physiologic model, with TCDD elimination as the main focus. The Dankovic model does not account for phenomena such as TCDD absorption, distribution, binding to liver receptors, enzyme induction, and synthesis of binding proteins, which occur on a much faster time scale (hours to days) than TCDD elimination (years in humans), nor for liver sequestration of TCDD (11,19). This reflects on the assumption of an-equilibrium between TCDD in lipid fraction of blood, liver, and adipose tissue. The model assumes that on a long-term basis the proportion of body mass represented by the adipose tissue becomes the major source of variation in TCDD kinetics across individuals and within individuals over time, given that adipose tissue in humans displays a much larger variation than liver volume.
The model, shown in Figure 1, is based on the assumption of a dynamic equilibrium of TCDD concentration between various body lipid compartments that form the total distribution volume (TLV= total lipid volume). The elimination ofTCDD due to liver degradation is assumed proportional to the total amount present in the liver with proportionality factor kf Moreover, a daily TCDD intake (pg/kg/day) proportional to body weight is assumed.
The distribution space of TCDD is partitioned into three subcompartments: adipose tissue, liver, and other tissues.
The actual distribution volumes (LV, grams) of TCDD are lipids, and are calculated for the above tissue compartments according to the International Commission on Radiological Protection (13), as 80, 6.9, and 2.2% of their respective  volume/weight (IV, grams). In particular: LVVother =0.022 Vother.
The MPTK model of TCDD is therefore described by the following linear, timevarying system with first-order dynamics: where the time dependency of liver lipid volume, total lipid volume, body weight, and daily intake has been indicated explicitly. In Equation 1, X(t) represents the total body TCDD in picograms (pg), with initial condition at time to given by Equation 2, which is calculated from the first measured lipid-adjusted serum concentration (ladj(to) (ppt)). Using the first data point as the initial condition makes it possible to disregard the exposure history before time to which then does not influence the TCDD dynamics after to. This approach was necessary because the exposure history for the RH individuals before to was not available. Equation 3 represents the prediction at time t of the lipidadjusted TCDD concentration. In the original model formulation (10), daily intake of TCDD was characterized by a parameter input describing background exposure per kilogram body weight. When applying the model to the NIOSH cohort, individual occupational exposure to TCDD is characterized as an additional daily intake per kilogram body weight proportional to the exposure time curve derived from the individual work history. In particular, intake(t) = input for the RH cohort, where input is a constant parameter, and intake(t) = input+ exposure up(t) [4] Environmental Health Perspectives * Vol 106, Supplement 2 * April 1998 for the NIOSH cohort, where u.,(t) is the exposure function and exposure is the unknown occupational exposure level (pg/kg/day), which is assumed to be unique for all exposed jobs. In the NIOSH cohort, the exposure function is represented by a list of time instants {Ito,t1,. . .. tJ and by a list of weights {we1,... ., we,j.
The weights wej take the values of 1 or 0 whether TCDD exposure has taken place in the time interval [t,.4, t.] or not. In mathematical terms, the exposure function for the i-th individual is defined as: uep(t,i)=wej(i) for tjl<t<t, and ue,(t,i)=0 for t < t0and t> t [5] A Reformulation of the TCDD Kinetics Model A useful simplification of the above model arises from the reasonable assumption that the body height of an adult subject does not change appreciably over time. With this assumption it is possible to rewrite the model equations taking into account only the BMI. This simplification can be useful in cases where individual body weight and height are unknown and the anthropometric characteristics need to be assigned using population values. In particular, it is more convenient to fix prior distribution and time variation only to BMI rather than to both body weight and H.
The simplification that restricts the applicability of the model to adults is based on the normalization of all quantities with respect to body weight (lower case will be used for the corresponding acronyms). For example, the average grams of adipose tissue per kilogram body weight are calculated as vadipose(t) = 1000 (0.01264 BMI(t)-0.13305), where BMI(t) is expressed in kg/mi2. The determination of other normalized tissues and lipid volumes is straightforward.
This simplified model was used to analyze both the RH and the NIOSH cohort data, using different descriptions of daily TCDD intake as described previously.

Comparison with Other ModelingApproaches
In the following we analyze the relationship between the above MPTK model of TCDD and another statistical, black-box, model proposed in literature, in particular the first-order kinetic model adopted by Michalek et al. (3).
It can be first noted that the dynamic Equation 6 has the following explicit analytical solution: x(t) = x(to)e Ito + It e inke(s) [9] where g(t) = kf liver(t) + (dBMI(t)ldt) ftlv(t) BMI(t) [10] represents the time-varying fractional clearance rate. Equation 9 follows from Equation 6 because linearity of the system has been assumed (20), i.e., the fractional clearance ofTCDD, g(t) given by Equation [ [11] where 'C(t) is the TCDD concentration t years after exposure, C0 is the initial concentration, and X is a constant but unknown decay rate" (3). Michalek et al. (3) corrected the TCDD values for background levels by subtracting 4 ppt. Parameter estimates were obtained after log transforming the data, i.e., for the model Equation 11 one obtains log(q t)-4)=log(G6)-Xt. [12] The relationship between the two modeling approaches arises by ignoring daily TCDD intake (i.e., intake(t) = 0) in Equation 9, and by taking the measurement Equation 8 into account. This yields the following equation: log(ladj(t)tlv(t)) = log(Iadj(to)tlv(to)) [13] which is similar to Equation 12, if one considers the equivalence X=g, with g representing the average value of fractional clearance, i.e., by putting to= 0 = fotg(r)dr t [14] In the interindividual variability of TCDD clearance, a statistical model was used in Michalek et al. (3), based on a mixedeffects linear approach both without and with adjustment for covariates. For the unadjusted case the model was [15] where subcripts i and j represent the subject and the sampling time, respectively, and (p and PI) represent the fixed population effects, ( ¶q and S-j) the random effects.
The adjustment for covariates, xT, was performed using the following model in Michalek et al. (3) log(Ci(tij) -4) = Ip+XtA + I tij+ 2 + P3Xijtij + £ij. [16] In the MPTK model (Equations 6-8), individual clearance and its variations is given by the integral term on the righthand side in Equation 13, which depends on changes in BMI (Equation 10). Interindividual variations of TCDD clearance are therefore assigned a priori. However, it must be stressed that this MPTK model allows for changes over time in measured lipid-adjusted TCDD concentration (given by Equation 8) independently from effective changes in total quantity of TCDD; this is because of the time-varying total lipid distribution volume.
The full model includes also daily TCDD intake, which consists ofbackground . It follows that kf is a nonlinear parameter while input and exposure are linear parameters for the function describing whole-body TCDD kinetics. In particular, the combination of Equations 9 and 4 with the output Equation 8 yields, for fixed functions g(t) and uexp(t), a linear relationship between TCDD predictions and parameters exposure, input, and even x(to). It is therefore possible to predict individual TCDD concentrations from BMI(t) and work history by a linear model, once the parameter kf has been fixed, i.e., solutions were constrained to go through the first data point, whereas error around the initial condition could be modeled with quantile plots of residuals (NLME the mixed-effects model. third Ranch Hand data point.
Naive Pooled Data Approach. The model parameters kf and input were ini-;-) represents the prediction tially estimated by NLWLS according to riven personal data history the NPD approach. This approach assumes Pi (which indudes BMI time that the estimated model parameters kf and sampling times, and work input are the same for all individuals. describes the initial TCDD Interindividual variability of TCDD clearnd ym, y2, and y3 are regresance is therefore attributed to changes in that can be obtained from BMI alone. To determine the effects on parameter dons estimates, we chose three different weight-lODs ing schemes: uniform weighting, reciprocal iations 6 to 8 were imple-of measurements, and reciprocal of squared nulated using the software measurements. From a statistical point of Numerical integration was view, these weighting schemes are optimal ig a fourth-order Runge-for the following assumptions on measurewith adaptive stepsize con-ment noise: constant variance; Poisson dis-I data sets, the time course tribution of measurements; and constant individual was assumed to coefficient of variation, respectively. In tween the sampling times, addition to the above weighting schemes, asures were only 5 years we considered also the NLWLS problem 4JIOSH cohort, a model of after log transformation of the data (and of BMI was implemented as the model output) using uniform weighting r. Simulations were carried Nonlinear Mixed-Effects Approach. of TCDD kinetics, we considered a nonlinear mixed-effects (NLME) model. We obtained simultaneous estimates of the population parameters kf and input as well as of the variance of the random effects associated with kf and with the measured TCDD concentrations (including the first measure which was used as initial condition for the simulations). The parameter estimation was carried out only with log-transformed data.
The approach was based on linearization of the model predictions (model output) for propagating the variability of the random effects as considered in (23). In particular, we assumed the following model for the fractional clearance parameter k yi)=k+ ekf (i), ekf (i)~NOc,f, [8 where i represents the i-th subject, k is the population mean and alf is the unknown variance of the random effect. For measurement noise, the following model was assumed for the log-transformed data Zi,j= log[ladj(kf (i), input, t()] + Ei(t) ££8~&, [19] where zt, represents the j-th observation of the concentration's logarithm for subject i. The fixed effects parameters kfand input, and the variance of the random effects a if and a2 were obtained by maximum likelihood estimation. The approximate covariance matrix of parameter estimates was calculated from the inverse of the Hessian matrix of the optimal cost function. Precision of parameter estimates was expressed either in terms of their standard deviation, computed as the square root of the corresponding diagonal element of the covariance matrix, or in terms of percent coefficient of variations, defined as 100 times the standard deviation divided by the parameter value.

Time Course
Availability of the time course of BMI is necessary to solve the MPTK model of TCDD. In the RH cohort, anthropometric measurements were available at intervals of approximately 5 years, such that linear interpolation of BMI measurements was acceptable for reconstructing individual time courses. On the contrary, the two values of BMI available for the NIOSH subcohort were in several cases measured decades apart. Therefore it was necessary to derive a model of BMI over time as a function of age. Data from the First National Health and Nutrition Examination Survey Epidemiologic Follow-up Study (24) provided preliminary information on mean and variance of the 10-year change in BMI by 10-year groups of age at baseline. The youngest group is the 25 to 34 years of age at baseline. To cover the spectrum of age at baseline induded in the NIOSH cohort, we included information on the additional decade 15 to 24 years of age at baseline. We obtained data on this age group from the Bogalusa Heart Study Group (25).
The overall time course of BMI variation resembles a descending staircase and the resulting BMI time course is a parabolalike piecewise linear function. However, this description did not appear to model accurately the BMI variations in the NIOSH subcohort. Therefore, we preferred to model time variations of BMI as a linear function of age using the subgroup of subjects with two BMI measurements. The adopted model was therefore dBMI =BMIt + PBMI [20] which yields the time course BMI(t) = BMI(t1) + aBMI (t2 -t2)

+ PBMI(t-tl )[211]
Estimates of aBMI and PBMI were obtained by fitting the squared difference between log-transformed BMI values observed at time t2 and values predicted according to Equation 21. The estimated values of OCBMI and PBMI were used in Equations 20 and 21 to solve Equation 6 for the computation of the exposure indices in the NIOSH cohort. In the subsample used to determine the occupational exposure parameter, we fixed the parameter aBMI to the above estimated value, and individualized parameter PIBMI to meet the second BMI measure.

Estimation ofOccupational Exposure
Individual occupational exposure levels of TCDD for a given work history, as available in the NIOSH cohort, were estimated from plasma concentration measurements taken one point in time. The estimation of the occupational exposure parameter was obtained via linear regression using Equation 17.
We estimated initially only parameter exposure, with both uniform weighting and after log transformation of the data, while ignoring the background exposure input. This was based on a preliminary multiple linear regression analysis in which the covariate y3 resulted the most important one for predicting TCDD concentrations. Despite a strong effect of measurement noise, we also evaluated the approximate individual exposure levels by the fraction ladjily3(t,pilkf), where dji represents the measured TCDD concentration in individual i.
The second application of Equation 17 was the simultaneous estimation of parameters exposure and inputwith the initial concentration ladj(to) fixed at 7 ppt. Again both uniform weighting and log transformation of the data were employed. The alternative approach of estimating all three parameters gave unrealistic estimates of ldj(to).

Computation ofExposure Indices
Computation of TCDD exposure indices from individual work histories was performed on the basis of simulated TCDD plasma concentration time courses as described in the previous sections. In particular, given an individual work history represented by Equation 5, and the time history of BMI(t) and its derivatives given by Equations 20 and 21, plasma TCDD concentrations were determined according to Equation 17 for fixed values of parameters 0 = [k.p input, exposure]. By dropping the explicit dependency on individual work history, BMI variations and parameters 0, we define, following Thomas (26), the general form of a cumulative exposure index computed at time Tfor the i-th subject, D( lT;n), as a weighted integral of the TCDD plasma concentration profile Di(T;n) = fto f(Tt,it)ldj(t, pi 0)dt, [22] where f(T-t, Ic) is a suitable weighting function parameterized by x, to represents the time of hire, and Tthe time at risk.
Typical choices of f(r,r) are a) the unweighted cumulative exposure with f(r,r) = 1, and b) the lagged cumulative exposure with f(er,i) = 1 if r. r, and flr,r)=0 if c<,r.
Sensitivities ofExposure Indices. Given that exposure indices are computed using population estimates of the kinetic parameters and of the occupational exposure, it is important to know how sensitive the computed indices are with respect to interindividual variations of the assigned model parameters. For this purpose, we observe that for a specific assigned model parameter Oj, the sensitivity of the exposure index can be computed as aDi(T;r) T aJad](t,pi I 0)dt. [23] aj .aladj(t,pi I 0) = y3(t, p, 0). [25] aexposure The computation of the sensitivity alad](t,pilO)/akf requires more complex calculations which were performed using the software PANSYM (21).

Parameter Estimation
Nonlinear Weighted Least-Squares. Parameter estimates are dependent on the weighting scheme and on log transformation of the data ( Table 4). The best model predictions and distribution of residuals (not shown, see below) were obtained with log transformation of the data, which also yields the smallest background input and the smallest value of the elimination parameter kf.
Nonlinear Mixed-Effects Approach. The parameter estimates obtained from logtransformed data with the nonlinear mixedeffects model are also reported in Table 4.
The random effect associated with the assignment of the first data point was taken into account in the analysis. Given the log transformation, a 2 is an estimate of the squared coefficient of variation of TCDD concentration measurements, which is then 25.9%. This value is likely to be overestimated and may reflect also interindividual variability of whose estimate, quantified by a k, resulted on the contrary to be very small. Nevertheless, the population mean estimate of kfis close to the values obtained with NLWLS and log transformation.
Model predictions and normal quantile plots of residuals are shown in Figure 2, separately for the second and third RH data points (although the fit was done simultaneously on all data points). No identifiable pattern was evident in plots of residuals versus fitted values for this model (not shown). These plots did not differ appreciably from those obtained under the NLWLS model and log transformation. The parameter estimates obtained with the NLME model and log transformation were used as reference values in the subsequent analyses.
Estimation ofOccupational Exposure. Estimates of the occupational exposure parameter were obtained by using Equation 17 with kf= 0.02199, corresponding to the value obtained with the nonlinear mixed-effects model and log-transformed data, and by fixing the initial condition of TCDD ladj(to) = 7 ppt (average level in an unexposed reference group of 79 workers) (17). The estimates obtained from the model ladj(t,pi kf) = ladj(to) Y1 (t,Pilkf) + exposure y3(t,pi) kf) (Equation 17) are shown in Table 5 for both uniform weighting and after log transformation of the data. In the same table we report the median and the 2.5 to 97.5 percentile interval of individual exposure levels. Although this is a rough estimate of individual TCDD exposure, because it is particularly sensible to measurement errors, it suggests that individual exposure levels are highly variable. Moreover, since the median value of this estimate is closer to the estimate obtained with log transformation of the data, it can be assumed once more that this latter estimation approach is preferable to ordinary least-squares. Results of simultaneous estimation of exposure and input parameters are reported in Table 6. Compared to Table 5, the estimate of exposure obtained with uniform weighting and with log transformation do not differ as widely.
The value of exposure used for computing exposure indices in the whole NIOSH cohort is the one reported in Table 6 obtained with log transformation. This choice was based on model predictions and on normal quantile plots of residuals (Figure 3), and it yielded the smallest parameter values, particularly regarding the input parameter. However, even the value of 0.45 for the input parameter would not be consistent with the assumed average population concentration of 7 ppt for unexposed subjects. In fact, assuming zero occupational exposure, the value of 0.45 would yield an average TCDD concentration of 10.3 ppt in the NIOSH subcohort instead of the postulated 7 ppt. We slightly adjusted this parameter to maintain an average concentration of 7 ppt. Therefore, for subsequent calculations we fixed background input at 0.293 (pg/kg/day).

Estimation ofPopuladon BMI Time Course
We obtained the following estimates: aBMI =-3.755 x 10-3 ± 0.9 X 10-3 (± SE) (kg/m2/year2), and PBMI = 0.26907 ± 0.04 (kg/m2/year). A comparison of the model predictions (Equations 20 and 21) with survey data regarding the rate of change of BMI over time and the corresponding variations from baseline is depicted in Figure 4.
Example of Calculation of Expsure Indices Figure 5 shows an example of application regarding the calculation of the time course of two exposure indices: serum dioxin concentration and its time integral (area under the curve). Each of these indices was calculated starting from age at hire and ending at age last observed. We also computed the sensitivities of these exposure indices to parameters: occupational exposure, background input, kf, and the assumed serum TCDD concentration at hire (initial condition). The sensitivities make it convenient to recalculate the exposure indices for different assumed values of the model parameters, without having to run additional model simulations. As an example, we report the recalculation of the cumulative exposure for a sizable variation (-30%) of kf and the corresponding approximation based on sensitivities. The approximation is fairly good for the nonlinear parameter kf. This approach yields exact results for deviations of any magnitude for the remaining linear parameters.

Discussion
We estimated occupational exposure to TCDD for members of the NIOSH cohort (14). Calculations were based on the kinetic model for TCDD proposed by Dankovic et al. (10). We first revised the model and worked out a simplified form based on the time course of BMI. We then carried out estimation of the parameters of this model (liver elimination constant kf (days-') and background input (pg/kg/day)) using data with repeated measures of serum TCDD taken over time (RH data and data on an unexposed reference group). Second, we  used the best estimates of the model parameters to estimate the occupational input rate in a subset of the NIOSH cohort for which single measures of serum TCDD were available. The occupational input rate thus estimated was then assumed to hold for all exposed jobs in the NIOSH cohort leading to a characterization of the time course of serum TCDD in individual cohort members, thus providing the basis for the calculation of exposure indices.
The elimination and is based on a minimal physiologic structure which includes the effects of variations in BMI. The model does not account for liver sequestration or binding of TCDD. Liver accumulation of TCDD may be of lesser relevance in the range of concentrations encountered in the RH group (Table 1). It may, however, be of importance for a sizable portion of the NIOSH data (Table 2). There is some evidence that human hepatocytes may be less sensitive than rat hepatocytes to the protein-inducing effect of TCDD (27). In the presence of a relevant liver sequestration of TCDD, its omission from the model might produce biased predictions, which, however, we did not observe.
The model displays a variation in TCDD kinetics through its dependency on changes in body mass and therefore in the lipid content of body compartments. This gray-box approach avoids the inclusion of covariates for statistical adjustment, for example the interplay of statistical covariates (body mass at a particular time, change in body mass over time and age) translates into a time-varying volume of distribution for TCDD. The MPTK model requires estimation of a smaller number of parameters than a statistical model with covariates and interaction terms. In addition, estimated parameters have a direct interpretation; and the availability of sensitivity indices to variations in model parameters constitutes a useful diagnostic tool. Prediction of TCDD serum profile over time and of derived quantities are straightforward even in presence of complicated exposure patterns such as work histories.

Parameter Estimation in the MPTK Model
Estimates obtained under the nonlinear mixed effects model (kf= 0.02199 days-', 95% CI=0.020, 0.024; input=0.1251 pg/kg/day, 95% CI=0.071, 0.179) were very consistent with the nonlinear weighted least-squares estimates obtained on log-transformed data and displayed the best behavior of model predictions and of model residuals. Nonetheless, the model fit was better for the second RH data point than for the third (Figure 2). With the NLME model we were also able to take into account the random effect associated with the assignment of the first data point. Estimates obtained under this model were used in subsequent calculations.
All estimates were obtained at the population level due to the sparse nature of the serum TCDD data used to fit the model. An analysis based on individual estimates (not shown) could only be performed for fixed values of the input parameter and it provided a similar estimate of kf albeit with a high dispersion.
Given the low values of TCDD measured in 1992 in many subjects and the reasonable precision in input estimates, the background subtraction approach, as considered by Michalek et al. (3), was not investigated. The reason for this was that by fixing an arbitrary lower bound for TCDD concentrations, we would have had to discard many subjects to maintain the final values positive, to be able to implement different weighting schemes, or to take the logarithm of the data.
Observations with nondetectable levels were excluded from our analysis. The exclusion affected nine observations in the reference group and three observations in the RH group, with the 1997 level less than 10 ppt. These observations would have provided additional information to estimate the background input parameter. However, imputing a value for these observations was complicated by the variability of the detection limit across measurements. This prevented the use of methods based on the assumption of a common detection limit (28). Some of the data used to estimate parameters in the MPTK model were the results of selection: for the exposed RH observations, the follow-up data on serum TCDD were available if the 1987 level was greater than 10 ppt. On the other hand a sample of RH veterans with 1987 serum TCDD less than 10 ppt was offered an additional measure in 1992. There were no selection criteria for the availability of serial measurements in the unexposed reference group. In the context of our modeling approach it was not possible to take into account these complex selection criteria from which the data arose. On the other hand, an analysis of the exposed Ranch Hand data based on a statistical model (3) did take the selection criteria into account following a data-conditioning approach. To assess potential biases in our analysis, we carried out a comparison of the predictions of the MPTK model and of the model in Michalek et al. (3). By taking hypothetical subjects with constant percent body fat over time and by setting the background input = 0, one can compare the apparent half-life of TCDD between the two models. Results in Table 7 show that within the conditions defined above the agreement between the two approaches is  have a clear impact on the apparent half-life ofTCDD. Figure 6 shows that higher values of BMI are associated with longer TCDD half-lives, whereas an increase/decrease in BMI in the same individual is accompanied by a decrease/increase in the apparent halflife of TCDD. This last effect, which is a distinctive feature of the model by Dankovic et al. (10) is due to a change in volumes of distribution. The curves in Figure 6 also represent system impulse responses following Equation 9, with x(to) = 100 and intake= 0.
The nonlinear mixed-effects approach used here may provide optimistic estimates of parameter variance. Other computationally intensive approaches (29) may provide a more realistic assessment of parameter variability at the population level.

Estimaon ofOccupation Exposure to TCDD
The estimation of the occupational intake rate was conducted by applying the MPTK model with kf = 0.02199 to the NIOSH subcohort of253 workers. We showed theoretically that there is a linear relation between TCDD serum concentration, occu-"pational exposure, background input, and initial TCDD serum concentration. On the other hand, the relation to kfis nonlinear. Estimates of the occupational exposure parameter were sensitive to data transformation. We selected the estimates based on logtransformed data, which yielded an occupational exposure rate of 232.7 pg/kg/day (95% CI = 192, 273). This choice was accompanied by the value of the input estimate, which was the closest value compatible with the observed average concentration of 7 ppt in absence of occupational exposure, and was supported by model predictions and residual plots. The estimate of background input of 0.45 pglkg/day was further adjusted to maintain a prediction of 7 ppt (input = 0.293 pg/kg/day). The need for this adjustment may indicate that parameter estimation in the NIOSH subcohort might benefit from the availability of an additional serial measurement of serum TCDD. While estimating occupational exposure in the NIOSH subcohort we had to assume that the occupational exposure intake of TCDD was identical across exposed jobs, given that a job-exposure matrix was not available. This has the effect of introducing a nondifferential misclassification of exposure, because of the absence of a relation with the outcome (disease) status. Although this has been traditionally associated with the introduction of a bias towards the null in the risk estimates (30), we feel that the direction of the bias is actually unknown, given that predicted serum TCDD is a continuous function of several variables and given the multivariate structure of the risk estimation models in which the TCDD exposure indexes eventually will be used (31).
The model fit to the NIOSH data showed a higher dispersion than observed in the RH data. This may be due to a combination of the following factors: a higher exposure level in the NIOSH cohort than in the RH group, differences in populations with possible effects on TCDD kinetics, the availability of a single TCDD measurement, and the assumption of a unique exposure level for all exposed jobs, as discussed above.
To account for BMI changes over time in the NIOSH subcohort, we did not rely on simple linear interpolation between the two data points, as we did with the RH data, for which the measurements were taken at relatively short time intervals. In fact, in the NIOSH subsample there is a large time difference between the first (at hire) and the second (several decades later) BMI measures. We selected a model structure for BMI changes over time that was compatible with BMI changes observed in survey data (24). We then estimated the parameters of this model using the data from the NIOSH subcohort. Given the important effects of BMI changes over serum TCDD kinetics time, we believe that the ad hoc model of BMI change should be more reliable than the use of general population survey data.

Computation ofExposure Indices
Finally, we calculated the time course of serum TCDD and of its area under the curve (cumulative dose) for individual members of the NIOSH cohort. This step was carried out with fixed values of the occupational exposure, background input, kf parameters, and of the assumed TCDD concentration at hire. In addition, this step requires knowledge of BMI at hire and of the complete work history. Each tabulated time of serum TCDD and of the area under the curve is also accompanied by sensitivity coefficients to occupational exposure, background input, kf, and the assumed TCDD concentration at hire. These sensitivity coefficients can be used to obtain alternative values for the exposure indices for different values of the parameters, without having to rerun the kinetic model simulations. The exposure indices thus recalculated are precise for sizable variations of kf (within 30%); they are, however, exact for deviations of any magnitude for the remaining parameters. We believe that the information on parameter sensitivities is very valuable in the dose-response analyses since it makes it convenient to build exposure indices for different values of the model parameters and then to evaluate the robustness of the risk estimates.