Population toxicokinetics of benzene.

In assessing the distribution and metabolism of toxic compounds in the body, measurements are not always feasible for ethical or technical reasons. Computer modeling offers a reasonable alternative, but the variability and complexity of biological systems pose unique challenges in model building and adjustment. Recent tools from population pharmacokinetics, Bayesian statistical inference, and physiological modeling can be brought together to solve these problems. As an example, we modeled the distribution and metabolism of benzene in humans. We derive statistical distributions for the parameters of a physiological model of benzene, on the basis of existing data. The model adequately fits both prior physiological information and experimental data. An estimate of the relationship between benzene exposure (up to 10 ppm) and fraction metabolized in the bone marrow is obtained and is shown to be linear for the subjects studied. Our median population estimate for the fraction of benzene metabolized, independent of exposure levels, is 52% (90% confidence interval, 47-67%). At levels approaching occupational inhalation exposure (continuous 1 ppm exposure), the estimated quantity metabolized in the bone marrow ranges from 2 to 40 mg/day.


Introduction
Benzene is carcinogenic in animals and humans via one or several of its metabolites (1). Thus, the fraction metabolized is likely to be a better measure of toxic exposure than benzene exposure itself. However, the extent of benzene intake metabolized at low exposure levels in humans is still unknown. For humans this fraction is difficult to measure directly, but can be estimated with a physiological toxicokinetic model (2)(3)(4)(5). These compartmental models ( Figure 1) allow the simulation of a variety of end points in specific organs, while providing the opportunity to use relevant prior information on physiological parameters, such as blood flows, organ volumes, etc. (6). Many of the model parameters, typically those controlling metabolism, are not known with precision. Proper statistical inference about the value of these parameters is therefore necessary, while conserving the strong prior information conferred by their physiological definition. Bayesian statistics also provide a natural way to merge a priori knowledge, gained by implementing a physiological model, with the experimental data. In addition, our interests reside in inference about benzene metabolism in humans, i.e., in a diverse population, rather than in any one individual studied in published experiments. We therefore designed a statistical model describing the relationships between individual and population physiological parameters to estimate population variability (7)(8)(9)(10). Linking a statistical model to a physiological mechanistic model may seem a daunting task, given the number of parameters involved. In fact, up to now, little consideration has been given to statistical issues when using such models (11); however, for such difficult problems, Bayesian numerical methods can provide solutions (10). We describe the application to our model of Markov-chain Monte Carlo simulation, which is a simple and powerful tool. As a result we report predictions and confidence bounds on various measures of benzene metabolism in a human population. We discuss how this information can improve our understanding of benzene toxicology.

Data and Models
The data consisted of the concentrations of benzene in exhaled air and venous blood, and phenol in urine, for three male volunteers exposed to benzene in an inhalation chamber during 4 hr (12). Data were collected during exposure and over the following 2 days. Two exposure levels were used: 1.7 ppm (5.2 pg/liter) and 10 ppm (30 pg/liter). Phenol concentrations were recorded only at the 10 ppm exposure level and were corrected for urine density. Two unreliable phenol measurements were discarded from this analysis (the corresponding urine density was low and density correction would not apply well). The small (barely detectable) contamination of the subjects' exhaled air prior to the 1.7-ppm exposure was ignored. This contamination cannot be a carryover from the previous 10-ppm exposure because 1 month separated the two experiments. In addition, the body weight of each individual was recorded (55, 73, and 90 kg, for subjects 1, 2, and 3, respectively), as well the minute volume for subject 2 (11 liters/min ± 10%). We used a physiologically based compartment i as a function of blood flow, pharmacokinetic (PBPK) model in which Fi, volume, Vi, arterial blood concentrathe human body is divided into five com-tion, CQ, and partition coefficient, Pi.
partments: poorly perfused tissues, well-These equations are linear, except for the perfused tissues, fat, bone marrow, and bone marrow and liver compartments, in liver ( Figure 1). These compartments are which a Michaelis-Menten term describes assumed to be homogeneous and distributhe metabolic clearance of benzene. For tion limited by blood flow. Pulmonary these compartments: exchanges are modeled by assuming instantaneous equilibrium between alveolar dC. FC (C C) V Cj air, venous blood, and arterial blood.  (12). The distribution given was used for the population extrapolations. dThe density of fat tissues is assumed to be 0.92. "This parameter was set at each Monte Carlo iteration so that the sum of the organ masses (plus skeleton, 17% of LM) matched the imposed or sampled body weight. eliminated from the liver and bone marrow to urine by a first-order process (with constant K.). Phenol itself represents a fraction, fp, of the urinary metabolites. Computing the concentration of phenol in urine requires defining the urine formation rate, F,. A list of the parameters is given in Table   1. The differentials were solved numerically using our own software, MCSim. The model allows us to compute, for given parameter values and initial exposure conditions, various quantities relevant for our purpose: concentration of benzene in blood or exhaled air, concentration of phenol in urine (during baseline conditions and benzene exposure), and quantity of benzene metabolized in a given period of time in the liver or bone marrow. The statistical model describing uncertainties and variabilities was constructed using a hierarchical population approach (9,13), as illustrated in Figure 2. It has two major components: the individual level and the population level. At the individual's level, for each of three subjects, exhaled air and blood benzene concentrations or urinary phenol concentrations (y) were measured experimentally. The expected values of these concentrations are a function (f) of exposure level (E), time (t), a set of physiological parameters of unknown values (0), and a set of measured, covariate parameters ((p). E, t, 0, and (p are subject specific. The function fis the nonlinear E.IFK.9. X2 Figure 2. Graph of the statistical model describing the dependence relationships between several groups of variables. P, prior distributions; p, mean population parameters; Y2, variances of the parameters in the population; E, benzene exposure concentrations; t, experimental sampling times; 0, unknown physiological parameters; (p, measured physiological parameters; f toxicokinetic model; y, measured benzene concentrations in blood or exhaled air; c&, variance of the experimental measurements.
Environmental Health Perspectives * Vol 104, Supplement 6 * December 1996 physiological model, described above. The concentrations actually observed are also affected by measurement errors, which are assumed, as usual, to be independent and log-normally distributed, with a mean of zero and a variance a2 (on the log scale). The variance vector a2 has three components, 6 2 for the measurements in blood, {2 for those in exhaled air, and (53 for those in urine, because these measurements have different experimental protocols and therefore are likely to have different precisions. Given its relatively large imprecision, the minute volume for subject 2 was considered as a data point, with a fixed variance on the log scale (the value chosen corresponds to a 10% coefficient of variation Three types of nodes are featured in Figure 2: a) square nodes represent variables for which the values are known by observation, such as y or (p; were fixed by the experimenters, such as E and t; or were fixed by ourselves, such as the prior on p and X2; b) circle nodes represent unknown variables, such as 0, 2, i, or X2; c) following the notation of Thomas et al. (14), the triangle represents the deterministic physiological model f. An arrow between two nodes indicates a direct statistical dependence between the variables of those nodes.

A Priori Parameter Distibutions
To take into account known physiological dependencies between the toxicokinetic model parameters (e.g., between organ volumes and body weight, or alveolar ventilation rate and cardiac output) several of them were linked to the lean body mass or other parameter values via scaling functions (15)(16)(17)(18) (Table 1). For example, volumes are input as fractions of the lean body weight, and the maximum rate of metabolism in liver as a power function of lean body weight. Note that cardiac output is computed as the sum of the organ flows. The scaling coefficients were the actual parameters used in input. At the population's level, we assumed that each component of the 0 parameter set is distributed log-normally, with population averages p and variances X2 (in log scale). We have some a priori knowledge of p and X2, at least in the form of standard values. We assigned a priori truncated normal distributions (with parameters M and S in log scale) to the population means p, and inverse gamma distributions (with parameter 0) for the population variances 402 We defined prior value for the hyperparameters M, S, and Io on the basis of the literature. The choice of values for these parameters and the bounds for truncation are summarized in Table 1. Truncations correspond for most parameters to 3 SDs to be subtracted or added to the mean. They are all biologically plausible. The exceptions are the two fractions V-b,bm and fp, for which such truncation would be meaningless, and the body weight. In setting uncertainties (but not for variabilities), we tried to be conservative and set the prior variances higher rather than lower when there was ambiguity in the biological literature (for example, with the metabolic parameters). This weights down the prior variances to "let the data speak." For convenience we give in Table 1  The values used for organ masses, when expressed as fractions of lean body weight, are usually considered as reference values for 35-year-old males (19,20). Volumes in liters and masses in kilograms have the same values since a density of 1.0 is assumed for all tissues, except for fat (density 0.92). Both the uncertainty on p and the heterogeneity of the fraction volumes in the population are estimated to be of the order of 10 to 15% (coefficient of variation), depending on the tissue group. As a consequence of scaling, the organ volumes are constrained by definition and have to sum to the lean body weight, not including bones, for each individual. We compute the volume of the poorly perfused compartment by difference so that the constraint is automatically satisfied at each iteration.
The geometric means of the perfusion rates per unit mass for the different compartments were set to usually accepted reference values (19,20). The mean ventilation over perfusion ratio (VPR) was set at 1.6 (21), since the subjects were allowed some activity after exposure. Exp(S) and exp(Xo) were set at 1.1 for the flows to various tissues, and 1.3 for VPR. This corresponds approximately to 10 and 30% variability, respectively. Using perfusion rates accounts for the covariance between total organ perfusion and organ weight.
The geometric means used for benzene blood/air and tissue over blood partition coefficient are taken from ranges and values previously published (22)(23)(24)(25). Partition coefficients may vary by a factor of 2 with hematocrit or depending on fasting, for example (26). Therefore, exp(S) was set at 1.5 for all partition coefficients. Exp(1o) was set to 1.3.
Prior estimates for the scaling coefficients for the population's maximum rates of metabolism in the liver and bone marrow, Vmaxi, Vmaxbm, and for the ratios of the maximum rates to the Michaelis-Menten coefficients, Kml, Kmbm, were taken from the literature (22,23,25). A large uncertainty is still associated with these numbers, and we chose a value of 5 for exp(S), except for Vmaxbm, which is the ratio of marrow to liver metabolism. We set exp(10) at 2. Thus, we believe these parameters to vary in the population of Pekari et al. subjects by about a factor of 2, but we are uncertain by a factor of 5 as to their population means. It would be difficult to express this sort of uncertainty without an explicit hierarchical model.
Prior variance for the average urine flow, Fu, was set on the basis of Rowland and Tozer (27) half-life of phenol, K, is 4.5 hr (28), so we set the excretion rate to 0.003/min, with a large uncertainty and a 50% CV for variability. The average amount of phenol in urine of subjects unexposed to benzene is 6 mg/liter and ranges from about 1 mg/liter to 40 mg/liter (12): given that at equilibrium Kf = (phenol concentration) x FPl/4, it can be deduced that the mean of Kf should be approximately 0.007 mg/min; we set uncertainty and variability to the plausible values used for K,.
At the individual level, we had no prior information for most of the parameters (except for the measured covariates), so information about the distribution of an individual's 0 parameter values is given by the experimental data and by the population parameters, p and X2.

Statistical Computations
A Bayesian analysis allowed us to combine two forms of information: "prior knowledge" from the scientific literature, and "data" from Pekari's experiments, in the context of the physiological compartmental model. Neither source of information is complete. If prior knowledge were sufficient, the experiments would not have had to be done, but Pekari's data alone are insufficient to pin down the parameters to reasonable values. Our goal was to fit the data using scientifically plausible parameter values.
The second interesting feature of the Bayesian approach is that it produces a Environmental Health Perspectives -Vol 104, Supplement 6 * December 1996 posterior distribution for the parameters, rather than a mere point estimate. Thus, the analysis outputs distributions of parameter values that are consistent with both the data and the prior information. Our statistical analysis yields distributional estimates (posterior distributions) of the parameters for each subject and for the population.
Current standard practice in Bayesian statistics is to summarize a complicated high-dimensional posterior distribution by random draws of the vector of parameters, in this case, from the distribution P(0, p, 2, &YIdata prior). The simulations can then be used to compute posterior distributions of estimands of interest, including individual parameters, and also derived quantities such as the proportion of benzene metabolized under specified conditions. Because 0 has many components, we use a combination of Gibbs sampling and Metropolis-Hasting Monte Carlo sampling to perform a random walk through the posterior distribution. These samplings are iterative procedures particularly convenient in the case of hierarchical models. They belong to a class of Markov-chain Monte Carlo techniques that has recently received much interest (8,10,(29)(30)(31). Five independent Monte Carlo runs were performed.
Convergence was monitored using the method of Gelman and Rubin (32).

Population Extapolation
The distributions of the fraction of benzene metabolized in the bone marrow or liver at various continuous exposure levels to benzene (0.001 ppm to 10 ppm) were obtained by simulation over 3 weeks. We verified that equilibrium was always reached in those conditions (the amounts metabolized over the last day differed at most by half a percent from the previous day). To compute the fraction metabolized over the last day, the amount metabolized was divided by the amount inhaled on the same day (i.e., the product of the alveolar ventilation volume for a day by the benzene inhalation level). These simulations were performed for the population by sampling one random parameter vector from N(p,l) for each of the 5000 estimates of p and S. This accounts for parameter covariance since the 5 x 1000 individual and population parameter sets are random draws from their joint (multivariate) distribution, not just from the marginal distributions, as would be the case in simple Monte Carlo simulations. For these simulations, body mass was also sampled lognormally (with geometric mean and standard deviation given in Table 1) (20).

Model Fit
The use of the Markov chain simulations, which reached approximate convergence in about 20,000 iterations, has allowed us to obtain a very good fit to the data of Pekari et al., while maintaining scientifically plausible parameter values. Figure 3 shows the data values predicted for each individual versus their observed counterparts (all data values are concentrations). Predictions were made with the highest posterior parameter values. This iteration is not much better than any of the last 5000 and is quite representative of the set. For an optimal fit, all points would fall on the diagonal, but such an adjustment is not expected given the analytical measurement errors in the data. The deviations here are small and the fit seems reasonable. Figures 4 and 5 give the simulated time profiles for subject 1, together with the data points. The experimental SDs estimates are 0.22 ± 0.016, 0.17 ± 0.013, and 0.22 ± 0.017 (in log space) for the venous blood, exhaled air, and urinary  In its present form, the model accounts for uncertainty and intersubject variability but not for intrasubject variability. The second exposure was performed 1 month after the first, and subject variability is the most likely explanation for the slight underor overestimations of the terminal slopes. A problem is that intrasubject variability confounds the dose effect. A better experimental design would have exposed the subjects twice to each benzene concentration level. The parameter values found for each subject should therefore be considered as approximate averages for that individual.

Posterior Distributions: Parameter Values
Among the results are the posterior distributions of all parameter values for individuals (whose precision is affected by measurement errors) and for the population (whose precision depends on population heterogeneity). The posterior distributions of the individuals' and population's parameters are summarized in Table 2 (the last 1000 iterations of the five runs were pooled, and the distributions are established with 5000 values).
The location of many parameters is noticeably different from the corresponding prior mean. However, the posterior distributions for most parameters are consistent with their prior distributions. The most important shifts are observed for the volume of fat, the partition coefficients and the metabolic parameters. The partition coefficients are quite close to those found by Watanabe et al. (25). However, the values of the scaling coefficients for the metabolic parameters Vmax and Km in the liver and bone marrow are different from what has been previously assumed or found. The mean of the scaling coefficient of the maximum rate of metabolism in the liver is about 60 times higher than the prior mean. The maximum rate in the bone marrow is expressed as a fraction, about 16%, of that in liver and therefore follows a similar trend. A tension between prior and data can even be noticed: the population mean for these parameters is lower that the individual values. This is because the prior constrains the population mean more than the individual values, which are more influenced by the data. The prior pulls down the estimated population average. The explanation of these findings is the following: as will be illustrated below, benzene metabolism is essentially linear in the dose range studied. The estimates of Vma and Km are therefore driven up by the estimation process. Only the Vma KKm scaling coefficients retain a meaning as they correspond to the first order rate of metabolism in the organs. The posterior mode (the point of highest probability) of the mean population Vmax1Km ratio for the liver is 0.6/min, with a population average of 0.76/min ( Figure 6). For the bone marrow the mode is 2.8/min and the population average 1.12/min, which indicates that metabolism is significant in this organ. VmaxlKm for the liver, min-1 Figure 6. Estimated distribution of the mean rate constant of benzene metabolism in liver for a human population. The histogram is based on 5000 Markov-chain Monte Carlo runs, and represents uncertainty in the mean value.
However, a very large uncertainty affects this value: the bone marrow is a "hidden" compartment (from the point of view of the data at hand), and little information is available about its metabolism in humans. Uncertainty (summarized by the SDs given in Table 2), however, is largely reduced (by comparison to the prior S) for most of the other parameters, indicating that the data brought important information about them. Interindividual variabilities are measured by the population standard deviations X, also given in Table 2. Those SDs correspond to factor 1.2 to 1.4 for the metabolic parameters. They are lower for the physiologic scaling coefficients. However, only three young white subjects were observed, and wider variations would certainly be found when observing a larger population. We noted above that it is likely that intrasubject variability affects the parameters (such as partition coefficients or metabolic parameters), but we did not attempt to estimate, as the data seemed insufficient for that purpose.

Posterior Distributions: Fraction of Benzene Metabolized
The relationship between inhalation exposure level and fraction of benzene metabolized per unit time in the bone marrow, after 3 weeks of continuous exposure, is presented in Figure 7 is obtained in the liver. This relationship is linear over the data range and beyond. However, the curve above 10-ppm exposure levels is an extrapolation into a region in which saturation could occur, and interpretation should be careful in this respect. Watanabe et al. (25) found that saturation could occur beyond 100 ppm. They also found slight nonlinear effects (S-shaped curving) in the relationship, which we do not see here. The data of Pekari et al. is of far better quality than those used by Watanabe et al., and we put more confidence in the present results. However, nonlinearity may still be found in other populations or if more subjects were to be studied. Further, linearity in primary metabolism does not imply linearity in the subsequent enzymatic transformations of the metabolites. Recent work by McDonald et al. (33) also support the notion that benzene metabolism is linear at low levels: see also the review by Smith (34). It is reassuring that, regardless of slight differences in the shape of the relationship, the quantities of metabolites predicted by Watanabe et al. and by us are close (geometric mean 10 mg per day with SD corresponding to a factor 1.6, at 1 ppm exposure). The population distribution of the total fraction metabolized is spread over a relatively narrow range. At low exposure (0.001 ppm), the mean of 5000 Monte Carlo estimates of this fraction is 57% (SD 6%) and the 5th and 95th percentiles are 47 and 67%, respectively ( Figure 8). The same results are obtained at high exposure (up to 10 ppm) because the relationship between amount metabolized and inhalation exposure is linear. It would be interesting to check what fraction of this variability can be explained by differences in P4502E1 activity, as measured, for example, by the chlorzoxazone assay (35,36). These results are indeed conditioned by the use of a small data set, with only three subjects of similar characteristics. While uncertainty could be reduced by additional analyses, population variability, which in this study is approximately as large as uncertainty about individual subjects, could increase when more subjects are included. This type of analysis requires a population pharmacokinetic approach and is more sophisticated than simple Monte Carlo simulations for uncertainty assessment. The results of the population approach are far more reliable because they rely as heavily on human data as on a priori physiological information. The method is of general applicability and provides a basis for statistically valid inference with physiological models.