Drinking Water and Health, Volume 8: Pharmacokinetics in Risk Assessment (1987)

Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

PART IV Uncertainties: Integration with Risk Assessment and Resou roes

Dealing with Uncertainty in Pharmacokinetic Models Using SIMUSOLV Gary E. Blau and W. Brock Neely INTRODUCTION In the course of investigating the behavior of a chemical, the toxicologist normally doses a mammalian species and follows the movement and dis- tribution of the chemical over a period of time. Having obtained this type of information, it is desirable to have a physiological explanation of the observed results. The usual method of arriving at such an explanation lies in characterizing the concentration-time data by the terms of a mathe- matical model. If the model is based on known mechanisms, it is possible to predict the behavior of the chemical in other species. The ultimate purpose of these types of studies is to use the results to assess the impact on other species, including man, to help quantify various risk assessment scenarios. The building of such models is an iterative process requiring a sophis- ticated mathematical analysis. This paper will discuss the art of model building and demonstrate how the recently developed software package SIMUSOLV (trademark of the Dow Chemical Company) can help the toxicologist analyze pharmacokinetic data in a relatively easy fashion. MODEL BU I LDI NO The first step in the building process is to define the problem to be solved. This can range from simply estimating the bioavailability of a 185

186 GARY E. BLAU AND W. BROCK NEELY chemical for a specific species to the elucidation of the fate and distribution in a variety of mammals, including man. Once the problem has been stated, the next step is to postulate several physiologically meaningful mechanisms to describe the set of data that has been collected. The third step is to use the data to discriminate between candidate models. To help make this choice the maximum likelihood principle will be introduced. Bayesian methods can also be used (Reilly and Blau, 19741. The utility of this later approach, however, can best be demonstrated for parameter estimation rather than model discrimination. Once the best model has been identified, classical statistical procedures can be used to measure the ad- equacy and to help select additional models. Quite frequently the existing data set does not contain sufficient infor- mation to discriminate rival models. In this case additional experiments must be designed and carried out. The type of new information required is determined by using the models to predict the behavior of the system under a variety of different conditions and to find that point where the predicted behavior is the most different. For example, suppose the problem is to determine whether a chemical follows a monophasic or a biphasic clearance. The key data points to discriminate these models would be at the latest time points that are still detectable by the analytical method being used. In this case the answer of where to look for additional infor- mation is obvious, but normally, this is not true. It is a big mistake to try and build models after the data have been collected. Model postulation and analysis should be part of the experimental process. SIMUSOLV has been developed to bring this sophisticated approach to the laboratory so that the experimentalist can quickly examine the data as soon as they are collected. Such instant feedback provides a greater breadth of knowledge with which to plan future experiments. After the best model has been identified and validated by a residual analysis and/or goodness of fit test, the next step is to determine the uncertainty in the parameter estimates. SIMUSOLV produces linear and nonlinear confidence regions for the parameters. In a similar fashion to the model discrimination discussed above, there may be insufficient data to produce acceptable confidence regions. If this happens, additional ex- periments must be designed to minimize the region of uncertainty of the parameters (Reilly and Blau, 19741. This entire description of the model- building scenario is summarized in Figure 1. Having gone through this extensive model-building process, the final task is the obvious one of using the model to solve the problem under consideration. Since model building is based on data, a word of caution is necessary. Data have errors associated with them. This uncertainty is an inherent part of both the model selected and the parameters estimated. Any use of the model must reflect this uncertainty in terms of confidence

PK Models Using SIMUSOLV 187 MODELING PROBLEM rPRIOR FACTS l | AND THEORIES ~ 1-' 1 NO /SUITABLE \ MODEL SELECTED,/ \/ YES ESTIMATE MODEL PARAMETERS < MODEL YES USE MODES PARAMETERS - FIGURE 1 Sequential model-building procedure. POSTULATE MODELS EXPERIMENTAL | I DATA l , _ COMPARED | DESIGN XPTS. TO MODELS I I DISCRIMINATE MODELS I _J NO EXPERIMENTAL DATA \ No al_ _ DESIGN EXPTS. TO REDUCE UNCERTAINTY

88 GARY E. BLAU AND W. BROCK NEELY limits in the values predicted by the model, as well as confidence regions around the estimated parameters. Finally, a model is developed to answer a specific question. There is always a great tendency to use the model for something other than its intended purpose. This is fraught with many dangers and should be avoided if at all possible. In the remainder of this paper reduction to practice of these philosophical concepts will be accomplished with SIMUSOLV. UNCERTAINTIES AND ERROR ANALYSIS IN MODEL BUILDING To illustrate the concept of model building consider the case of two models (M~ and M2) characterized by the following equations: For M1: hi = f ~ Ski ~k2,xi) + c For M2: Yi = f2(k~,k2,k3,xi) + c (1) Each model is a function of two or more parameters, i.e., rate constants for clearance, metabolism, absorption, etc., and a single independent variable (xi), usually time. In both cases Hi is the dependent variable, such as concentration in the blood, and ci is the experimental error in the ith measurement. The first problem in uncertainty is to examine the two models and apply a set of criteria (to be discussed below) to make a decision as to which one is best. In other words, find the most suitable of the models described in Equation 1 which best describes the data. This is done by comparing how well the models fit the data. The comparison is done when the models are examined at their individual best. The actual procedure is accomplished by measuring the probability of generating the observed concentration- time data by each model for a given set of parameters. The method of maximum likelihood accepts it as obvious that the values of the parameters k~,k2 for Me and k~,k2,k3 for M2 which maximize this probability are the best. These parameters are called maximum likelihood estimators denoted k~*,k2* form and k~*,k2*,k3* for M2. The function which calculates the probabilities for any set of parameters is called the likelihood function, denoted as L~(k~,k2) for Me and L2(k~,k2,k3) for M2. To distinguish one model from another, the ratio of likelihood functions is evaluated by using the maximum likelihood estimators. A ratio of 10 is ordinarily taken as showing a difference in plausibility, whereas 100 denotes a strong pref- erence for one model over the other (Reilly and Blau, 19741. For example, if L~(k~*,k2*11L2(k~*,k2*,k3*) for Models 1 and 2 were 100, there would be a strong preference for Me being better at explaining the data than M2. The above heuristics assume that the number of parameters in each model

PK Models Using SIMUSOLV 189 is the same. If the number is different, then the likelihood ratios must be somewhat higher than 100 to discriminate between them. The above illustration can be generalized to M models as follows. For M candidate models of the form: Hi = fn~kn,ti) + e~ (n = 1,2,. . . M), (2) where k becomes a dimensional vector of parameters for the nth model selected. For each of the M models, SIMUSOLV can be used to calculate the maximum likelihood estimates kn*. From these estimates likelihood ratios can be determined for all pairs of models. The heuristics can be applied as described above in the two-by-two comparison test, and the most suitable model can be selected. An example of this technique is described in the section "Applications." The second problem in uncertainty is obtaining the best estimate of the parameters for the model selected. When the method of maximum like- lihood is used, this task is embedded within the discrimination problem as described above. In other words, by comparing the models at their optimum the best set of values is automatically generated. The determi- nation of these maximum likelihood estimates is a nontrivial task, except for the special case in which the models in Equation 1 are linear. Unfor- tunately, this situation does not occur with phenomenologically based models; consequently, it is necessary to apply estimation procedures (Bard, 19741. SIMUSOLV uses both direct and indirect algorithms for parameter estimation. In the direct method different values of the parameters are selected, and the likelihood function is estimated for each. The set which gives the lowest value of the likelihood function is deleted, and a new point is selected according to a procedure called the flexible polygon method (Bard, 19741. This process of adding and deleting points is continued until the maximum value of the likelihood function is obtained. Indirect methods not only use values of the likelihood function but also its curvature to direct the search for the maximum value. In SIMUSOLV, the generalized reduced gradient method is used (Lasdon et al., 19781. From an initial guess of the maximum, succeeding estimates are selected along a direction of improved function values. Determination of this di- rection comes from a knowledge of the gradient, which can be determined numerically or analytically. The former is usually preferred because the user does not need to supply any information. Numerical approximations to the gradient, however, frequently lead to less than adequate directions. This problem has been overcome in SIMUSOLV by solving both the sensitivity and the model equations to yield the gradients directly. This method, called the direct decoupled method, has reduced the computa- tional time over the conventional reduced gradient method by an order of

190 GARY E. BLAU AND W. BROCK NEELY magnitude (Dunker, 19841. In addition, the state-of-the-art numerical in- tegration techniques Lawrence Solver for Ordinary Differential Equations (LSODE) are being used to solve all the differential equations associated with the model (Hindmarsh, 19821. Despite the sophistications of the techniques included in SIMUSOLV, every effort has been made to make these capabilities invisible to the user. However, the user must be aware that the true maximum likelihood estimates require this level of technology to ensure efficient use of the computer resources and achievement of the proper solution. Once point estimates are obtained by maximizing the likelihood func- tion, their uncertainty must be quantified. SIMUSOLV uses several meth- ods to handle this quantification. The first method, called linearization, represents the region of uncer- tainty as an ellipse in parameter space (e.g., two-dimensional contour plots for the two parameters kit and k2 in Model 11. Here, the likelihood function is linearized by a Taylor's series expansion about the maximum likelihood estimates. These confidence regions are drawn based on a statistical F distribution (Draper and Smith, 19811. A 95% confidence region means that if the experiment was repeated 100 times, 95 times the maximum likelihood estimates would fall within this region. The linear- ization method is attractive in that the confidence regions are easily gen- erated by the computer (Reilly, 19761. Unfortunately, linearizing the function distorts the shape of the confidence regions. Rather than ellipses, confi- dence regions may take on different shapes dependent on the information content of the experimental data. The second method available in SIMUSOLV to determine these con- fidence regions, called the nonlinear method, does not linearize the like- lihood functions but uses the familiar F statistic to determine the region. In Problem 2 in the section "Applications," it will be shown just how significant these departures from linearity can be in the shape of the regions. The determination of the confidence regions requires a large amount of computational effort because the likelihood function must be evaluated for all points in parameter space. It is also an approximate method because the F statistic, which is being used, is valid only for linear models. The third method in SIMUSOLV for representing the confidence re- gions, called the exact method, requires no assumptions. It is based on a subjective interpretation of statistics called Bayesian methods (Reilly and Blau, 1974~. By this method the degree of belief in an event is quantified versus the familiar objective or frequency of occurrence interpretation. The Bayesian approach applied to the quantification of parameter uncer- tainty is expressed by the following theorem: Privily) = Prick) PrLylk), where Prickly) is the posterior probability of the parameter k and is the

PK Models Using SIMUSOLV 19i true value given the set of data y, Prick) is the prior probability of the parameter k and is the true value before the data were collected, and Pr (ylk) is the likelihood of generating the data y for the set of parameters k. Although the above likelihood is different in interpretation, it gives the same values as the likelihood function discussed earlier. To apply the method, it is necessary to multiply the believed value of the parameters before the experiment is conducted by the value of the likelihood function calculated after the experiments are performed. These values must be normalized so that the total values of the posterior probabilities add to 1.0. If no prior information is available for the parameters before the experiment is conducted, a uniform, noninformative prior distribution is assumed by SIMUSOLV (Bard, 19741. By using this posterior probability, SIMUSOLV calculates regions of uncertainty by using frequency of oc- currence arguments independent of any assumptions in the form of the likelihood function. Once again, an example later in the text will illustrate the effect of the assumptions inherent in the last two methods. SIMUSOLV allows the user to construct confidence regions by either of these methods. The computational burden associated with the nonlinear and exact methods is prohibitive for large problems (where the number of parameters is greater than four). They must be used, however, if there is any suspicion about the inflation content of the data. The Bayesian method should be used exclusively if prior information on the parameters is available. Experimental Error Before concluding this section, the effect of experimental error on model development must be discussed. Because the method of maximum like- lihood is being used, it is necessary to accommodate the structure of the experimental error in the model-building process. To simplify the analysis, it is frequently assumed that the errors in different experiments are (1) independent of one another, (2) uncorrelated, (3) have constant var- iance, and (4) are normally distributed. If these assumptions are not rea- sonably valid it may have serious effects on the statistical analysis. To illustrate, suppose that blood plasma concentrations are measured by an analytical procedure in which the error in the measurement is a constant fraction of the quantity being measured. Let the model expressing the blood plasma concentrations as a function of rate parameters be as follows: Cj = f~k~,k2, · ·, ti) + ei, (3) where the kj's are the individual rate constants, and the error ci is a constant fraction of the amount being measured. The variance of the error in the observations will not be constant if the blood level extends over any

|92 GARY E. BLAU AND W. BROCK NEELY considerable range. On the other hand, if the assumption is made that the variance of the error is constant, the consequences both for parameter estimation and model building could be equally serious. For example, during a long clearance phase from the body, the model will fail to give proper weight to any tailing in the curve from secondary processes. This will be demonstrated more clearly in the next section ("Applications") when the data in Problem 1 are analyzed. In the past, much effort has been expended in getting models into convenient form for plotting on ordinary or special graph paper. A bio- logical example is the use of Lineweaver-Burk plots for determining kinetic constants in enzyme systems. In most cases this convenience is achieved at the cost of distorting the error structure. When curves are fitted by eye or by simple linear least squares, constant variance is almost invariably assumed. For example, one method of plotting blood clearance data is by using logarithmic graph paper. This assumes that the elimination from the blood follows a first-order process and can be expressed as: dCldt = —kt, (4) where C is the blood concentration, and k is the elimination rate constant, which on integration yields the following nonlinear equation: C = C0 expf - kt) * c, (5) where ~ is the error in measurement and is a multiplicative function of the concentration. By making the log transform of Equation 4, a linear equation results. In this case the error becomes an additive function of the log of the concentration, as shown in Equation 6: In C = —At + In C0 + In e. (6) Plotting the logarithm of the concentration C versus time yields a straight line. The straight line, however, is achieved by assuming that the scatter in the data is a constant fraction of the quantity being measured when, in fact, it might be that the scatter in the data is constant. If it is the latter, then the width of the error bars as the concentration becomes small will make the interpretation of a straight line from Equation 5 misleading. It is generally much better to write the model directly in the form of Equation 1 and use SIMUSOLV, which can accommodate the error struc- ture. It is more the rule than the exception that the error, while unpre- dictable, depends to some extent on the magnitude of the quantity being measured. That is, the absolute value of the error usually tends to be large when larger quantities are measured. In the SIMUSOLV program the error is described by the following model (Reilly et al., 19771: Variance ~ Hi) = Variance (e ~) = o2 Ef (k,xi) ~ A (7)

PK Nlodels Using SIMUSOLV where ~ and By are statistical parameters of the error model. The latter is called the heteroscedasticity parameter. Usually its value will be between 0, in which case the error variance is a constant, and 2, when the variance is proportional to the quantity being measured. Equation 7 is used to weight the measurements to allow for changes in error variance. Although the maximum likelihood approach is restrictive in requiring explicit assumptions concerning the form of the experimental error, any unknown parameters such as ~ and By appearing in the error model are estimated along with the pharmacokinetic parameters of the model. APPLICATIONS The SIMUSOLV computer program was used to analyze the three simulated problems presented below. This is a package that has been developed at the Dow Chemical Co. to aid the nonmathematician in build- ing and solving models containing systems of algebraic and ordinary differential equations. After the program is wntten, SIMUSOLV allows the investigator to match experimental data with different models. Based on the statistical discussion in the previous section, guidance is provided in deciding which model is "best." If the information content of the experimental data is inadequate to make such decisions, SIMUSOLV will help the user design additional experiments to make the final model se- lection. Finally, SIMUSOLV contains graphical routines to display sim- ulated results. Problem ~ Assume that an animal has been given a C~4-labeled chemical by means of a single oral dose. At periodic time intervals blood samples are collected and analyzed for total radioactivity. The data in Table 1 are the result of such a hypothetical expenment. TABLE 1 Simulated Sequential Blood Concen- tration-Time Data Following an Oral Dose of 45 mg Time (h) 0.5 1.0 1.5 2.0 3.0 4.0 6.0 8.0 10.0 Blood (mg/liter) 0.10 0.11 0.082 0.057 0.035 0.024 0.0091 0.0050 0.0036

i94 GARY E. BLAU AND W. BROCK NEELY Two possible mechanisms to explain the data in Table 1 are the fol- lowing: A. A simple two-compartment model as shown below: ml kl2 k2o > 1 C2 1 3 U v2 Model A where m1 is the mass of chemical in the gastrointestinal tract at time t, V2 is the volume of distribution for the bloodstream (200 ml), C2 is the concentration of chemical in the bloodstream, and u is the mass of chemical in the urine. Assuming a first-order absorption of the oral dose from the gastroin- testinal tract into the blood stream followed by a first-order elimination into the urine, the mass-balance equations describing Model A become: dmlldt = —kl2m1 dC21dt= kl2(ml/V2)—k20C2 duldt= k2oc2v2 (8) B. A more sophisticated model can be postulated with the scheme shown below: ml kl2 > k2o c2 1 > u v2 k23 ~ ~ k32 m3 V3 Model B

PK Models Using SIMUSOLV 195 TABLE 2 Calculations from a SIMUSOLV Anaylsis on the Two Models A and B Heteroscedasticity Likelihood Percent Variance Model Parameter Function Explained A 0 1.0 x 10~6 97.66 B 0.89 1.4 x 102° 99.56 where m3 is the mass of chemical in compartment V3 that is not subject to elimination. In this case the parent compound is capable of moving into a com- partment that is not subject to elimination, a so-called deep compartment. Such a movement has the effect of lowering the amount of chemical that can be analyzed in the blood. The slow rate of release (k32) back into the blood generates a long clearance time into the urine pool (k32 ~ k201. The following set of differential equations describes this model: dm,,ldt = —karma dC21dt = k~2 m~lV2 + k32m31V2 k23C2—k20C2 dm31dt= k23V2C2—k32m3 duldt= k20C2V2, (9) where k32 has units of liters per hour and represents the volumetric transfer rate coefficient for the movement of material from the deep compartment (V3) back into the plasma (V21. With a little imagination more models could be postulated to explain the data in Table 1. These two models suffice, however, for our purpose, which is to address the problems associated with discriminating models. By plotting the logarithm of the concentration versus time, a straight line can be drawn by eye through the points, as illustrated in Figure 2. Without any further consideration, it might be concluded that the simple elimination as described in Model A would explain the results. By using likelihood ratios, however, discrimination between the two models is possible. An example of the complete output from SIMUSOLV for Model A is included in Appendix A. The key statistics used for discrimination are shown in Table 2. The likelihood ratio of 14,000 indicates that there is a strong preference for Model B over Model A. This conclusion could have been missed if

196 GARY E. BLAU AND W. BROCK NEELY 0.1 o m 0.01 z z o F E z O 0.001 o Cal Hi\ - - \ \ \ _ 1 1 1 1 2 4 6 8 HOURS FIGURE 2 Log concentrations versus time for the data in Table 1. 10 reliance on the simple log plot shown in Figure 2 had been made. The actual comparison of the values predicted by the models using their max- imum likelihood estimates and the data are shown in Figure 3. Here, it becomes quite evident that Model A does not describe the blood concen- tration points at the later time intervals. Thus, SIMUSOLV quickly allows the investigator to model the data; and by examining the key statistics and plots, the investigator can begin the task of eliminating and targeting in on the mechanisms that need to be examined in greater detail. Problem 2 Problem 2 demonstrates the ability of SIMUSOLV to examine the confidence regions around the rate constants. Metzler (1968) presented data for drug absorption and excretion (Table 31. The problem was formulated in the manner described for Model A in Equation 8. Figure 4 shows the agreement between the model and the data by using the maximum likelihood estimators. The information content

PK Models Using SIMUSOLV 197 0.12 0.10 - z o ~ 0.07 E an O 005 a: in A 0.02 0.15 - ~ 0.12 - z o 0.09 of O 0.06 a: co 0.03 0.001 A_ 1 0.00 ~1 o.oo, L I 0.00 1.50 3.00 rat ·\ \. MODEL A 4~50 6.00 7.50 TIME (hours) 9.00 10.5 MODEL B - - 1.50 3.00 4.50 6.00 7.50 9.00 10.5 TIME (hours) FIGURE 3 Concentration-time data from Table 1 fit to Models A and B. associated with the data for generating the uptake rate constant is very minimal (one datum point at 1 h from Table 31. This would indicate a strong suspicion that the confidence region around kl is very large. By using the linearization approach for quantifying the confidence regions, the contour plot shown in Figure 5 is obtained. Reliance on this approach

|9~3 GARY E. BLAU AND W. BROCK NEELY TABLE 3 Sequential Blood Concentration- Time Data Following an Oral Dose Time (h) Concentration 2 4 7 24 48 55 72 70.2 81.9 76.4 74.1 50.3 32.4 20.5 18.3 Source: Metzler (1968). would lead to the conclusion that the confidence region around the two rate constants is excellent. The distortion of the region is readily apparent when the exact procedure is used to depict the confidence region (Figure 61. As suggested by the data, the exact confidence region demonstrates the large uncertainty associated with the uptake rate constant. Figure 7 indicates the result of taking data points at 0.5 and 0.75 h following administrations. 100 80 o it 60 Cal he o <) 40 o o m 20 o - - - - - - . - - 1 1 1 1 1 1 1 20 30 40 50 60 70 80 0 10 TIME (hours) FIGURE 4 Fit of the data in Table 3 to the model described by Equation 8.

PK Models Using SIMUSOLV 199 TABLE 4 Plasma and Urine Concentrations Resulting from an Oral Dose of 48.15 mg Plasma (mg/liters) Urine (mg) Time (h) Parent Metabolite Parent Metabolite 0.82 0.175 0.822 1.0 1.87 7.23 1.2 0.166 1.14 1.4 0.126 1.15 2.0 0.109 0.860 3.23 15.53 2.4 0.09 0.648 2.9 0.083 0.601 3.0 4.02 21.15 3.38 0.070 0.382 3.92 0.059 0.403 4.0 4.59 25.88 4.42 0.051 0.304 5.18 0.035 0.252 6.0 5.77 32.4 6.35 0.015 0.143 8.0 6.3 34.9 8.3 0.0081 0.0636 10 0.0047 0.0332 12 6.65 37.06 12.4 0.0026 0.0065 24 6.92 38.7 24.57 0.0009 0.0065 48 7.3 40.29 72 7.38 40.77 The confidence region has now been reduced, and the new value for k1 is more reliable. Thus, by using the graphical capabilities of SIMU- SOLV, the choice of experiments for improving the quality of the data is greatly improved. Problem 3 Nichols and Peck (1982) orally dosed an experimental animal. Plasma concentrations of the parent and the metabolite were measured up to 24 h after dosing. Urine collections were made over 72 h, and the cumulative excretion of parent and metabolite in the urine was determined. The data for this experiment are shown in Table 4. Several models were postulated to explain the data ranging from the simple (Equation 10) to the complex (Equation 13~.

200 GARY E. BLAU AND W. BROCK NEELY CONTOUR PLOT 0.0 kit 2.0 4.0 ,,'"o j,~35~0o$~"` L ~^,,~r 0 lo ,~/ 1 1 (9 ~ al 1 .o 1 1 1 0.021 0.024 0.027 0.030 0.015 0.018 k2 FIGURE 5 Linear confidence regions around kit and k2 for Model A using the data in Table 3. Model 1 m kl2 c2 V2 1 k24 C4 (5 parameters) k2o , U2 k40 ~ U4 (10) where m1 is the oral dose, kl2 is the absorption rate constant, k20 is the excretion constant for parent, k40 is the excretion constant for metabolite, C2 is the concentration of parent in plasma, C4 iS the concentration of metabolite, u2 is the cumulative amount of parent in urine, U4 iS the cumulative amount of metabolite in the urine, and V2 is the volume of distribution for plasma. Note that all rate constants are first order.

PK Models Using SIMUSOLV 201 CONTOUR PLOT 0.0 2.0 k1 * / ~~'~ 70.0 ~S~;~` 4.0 6.0 8.0 ' ° ~ I '1N 11 1 (9 1 1 0 ' 1 1 i (D ~ 0 ' ''' J ° 1 CD in A l 1 0.015 0.018 0.021 k2 0.024 0.027 0.030 FIGURE 6 Confidence region around k, and k2 for Model A using the data in Table 3. Note the large uncertainty in the confidence around the uptake rate constant.

202 GARY E. BLAU AND W. BROCK NEELY CONFIDENCE INTERVAL PLOT 0.0 k1 2.0 <~2," 0.015 0.018 0.021 0.024 k2 0.027 0.030 FIGURE 7 Confidence region around k' and k' for Model A using the data in Table 3 with additional data taken at 0.5 and 0.75 h. Model 2 ml kl2 k2o , C2 1 > u2 v2 1 k24 k40 V4 C4 1 > U4 (6 parameters) (1 1) where all the symbols are similar to those in Equation 10, with the ad- ditional parameter of V4, the volume of distribution of metabolite in the plasma.

PK Models Using SIMUSOLV 203 Model 3 ml kl2 c2 \k24 \ C4 k23 ~ ~ k32 m3 k2o - > U2 k40 ~ U2 (12) (7 parameters) where two additional parameters k23 and k32 that are similar to those in Model B (Equation 9) are included. Again, these two constants represent the flow of chemical into a deep compartment (V3) and the slow release back into the plasma, or compartment V2. Model 4 ml m3 , //~;, 1 kl2 k2o > c2 - > U2 c2 k24 V4 k40 U4 > (13) (8 parameters) This is the complete model in which the volume of distribution for the metabolite (V4) has been added to the model represented by Equation 12. SIMUSOLV was used to determine the maximum likelihood estimates shown in Table 5. Table 5 also shows the progressive improvement in the fit of the model in going from Model 1 to Model 3. By applying likelihood ratios in a pairwise fashion, no discrimination is possible be- tween Model 3 and Model 4. However, Model 2 is slightly better than Model l, and Model 3 is significantly preferred over Model 2. The lack of discrimination between Models 3 and 4 is also evident when the volume of distributions V2 and V4 are examined. In other words, there is no

204 GARY E. BLAU AND W. BROCK NEELY TABLE 5 Statistical Analysis Using SIMUSOLV on the Four Models Represented by Equations 10-13 and the Data in Table 5 Model 2 4 Maximum likelihooda 57.34 66.28 90.15 91.65 Parametersb kit 0.362 0.365 0.36 0.34 k24 6.85 5.08 17.37 17.15 k20 1.11 0.83 3.08 3.06 k40 1.37 2.68 2.96 3.14 k23 3.34 2.48 k32 0.069 0.047 v2 9.8 13.6 3.33 3.41 V4 4.19 3.19 Percent variation explained 98.79 98.77 98.76 98.75 aThe maximum likelihood function is expressed in natural logarithms. Thus, the ratio of the function is simply the difference. ball rate constants are in units of reciprocal hours, except for k32, which has units of liters per hour. difference between the volume of distribution for the parent and the volume of distribution selected for the metabolite. The big improvement in fitting the data came from the addition of k23 and k32 to Model 3. The movement of the parent compound into a deep compartment and the slow release back into the plasma allowed for the improved fit at the later time intervals. As can be seen from the plot in Figure 8(A), it is the later time intervals that are not adequately fitted by Model 1. Figure 8(B) represents the fit of Model 3 and is much improved for all time periods. CONCLUSION The purpose of this paper was to demonstrate areas of uncertainty which arise in building pharmacokinetic models and the technique required to deal with them. In discriminating rival models and quantifying uncertainty in the model parameters, the inherent error in the actual experimental data must be considered. SIMUSOLV, a computer software package developed at the Dow Chemical Co., provides a tool for the nonmathematician to work with to handle these problems. By using the latest statistical tech-

PK Models Using Sl MUSOLV 205 60 - ~ 45 a) ._ ~ 30 o E - ce o o 60 ce ca 45 Q ._ ~ 30 o ct ct ° 15 o · 1 / - . . (A) _ ~ r (B) . TIME (HOURS) FIGURE 8 A plot of the data in Table 5 fitted to Model 1 (Equation 10). Note the divergence of the lines from the data at the later time periods. The lower curve is the fit of the same data to Model 3 (Equation 12). The incorporation of the additional mechanism has caused a much better fit to the data at all time periods. niques and employing a graphics package, the user is aided in making the proper choice between models. In addition, the package also provides assistance in the area of experimental design. In other words, where the information in the data is insufficient to make a choice, the program will suggest areas in which additional experiments should be conducted. REFERENCES Bard, Y. 1974. Nonlinear Parameter Estimation. New York: Academic Press. Draper, N. R., and H. Smith. 1981. Applied Regression Analysis, 2nd ed. New York: Wiley-Interscience.

206 GARY E. BLAU AND W. BROCK NEELY Dunker, A. 1984. The decoupled-direct method for calculating sensitivity coefficients in chemical kinetics. J. Chem. Phys. 81:2085-2093. Hindmarsh, A. C. 1982. ODePack, A Systematized Collection of ODE Papers. Report no. UCRL-880077, August 1982. Livermore, Calif.: Lawrence Livermore National Labo- ratory. Lasdon, L. S., A. D. Waren, and M. Rather. 1978. GRGZ User's Guide. Cleveland: Department of Computer and Information Services, Cleveland State University. Metzler, C. 1968. Data presented at the Pharmaceutical Science Meeting, November 17- 20, Washington, D.C. Nichols, A. T., and C. C. Peck. 1982. VIII a non-tr~vial example. Pp. 24-37 in ELSNR- Extended Least Squares Nonlinear Regression Program Users Manual, Technical Report no. 9. Bethesda, Md.: Division of Clinical Pharmacology, Uniformed Services Univer- sity. Reilly, P. M. 1976. The numerical computation of posterior distributions in Bayesian statistical inference. J. R. Stat. Soc. Ser. C Appl. Stat. 21:201-209. Reilly, P. M., and G. E. Blau. 1974. The use of statistical methods to build mathematical models of chemical reacting systems. Can. J. Chem. Eng. 52:289. Reilly, P. M., R. Bajramovic, G. E. Blau, D. R. Branson, and M. W. Sauerhoff. 1977. Guidelines for the optimal design of experiments to estimate parameters in first order kinetic models. Can. J. Chem. Eng. 55:614. APPENDIX The following is an example of the statistical output for Model A and is typical for all SIMUSOLV runs. Standard Description Parameter Estimates Deviation Initial Final Objective Function 32.005 36.876 k,2 1.7820 1.5317 0.186 k20 0.7899 0.9693 6.39E—02 Time C, C2 Standardized Residual Observed Predicted % Error Residuala Plot 0.00E + 00 0.0000E + 00 0.0000E + 00 0.00 0.0000E + 00 0.500 1.0000E - 01 9.8685E - 02 1.31 1.315E - 03 * 1.000 0.1100 0.1067 3.03 3.336E - 03 ** 1.500 8.0000E- 02 8.7028E - 02 - 8.78 - 7.028E- 03 **** 2.000 5.7000E - 02 6.3420E - 02 - 11.44 - 6.520E - 03 **** 3.000 3.5000E-02 2.9080E-02 16.91 5.920E-03 **** 4.000 2.4000E—02 1.2109E—02 49.55 1.189E—02 ******* 6.000 9.0000E—03 1.8813E—03 79.10 7.119E—03 ****** 8.000 5.0000E - 03 2.7726E - 04 94.45 4.723E - 03 *** 10.000 3.6000E-03 4.0222E-05 98.88 3.560E-03 **

PK Models Using SIMUSOLV 207 Statistical Summary Maximized Weighted Residual Weighted Likelihood Sum of Residual Function Squares Sum Standard Error of Estimate Percentage Variation Explained Weighting Parameters 97.667 0.00 36.88 3.668E - 04 Correlation Matrix kl2 1.000 0.1293 kl2 k2o Variance-Covariance Matrix kl2 k2o 1 .000 kl2 3.4559E—02 1.5384E - 03 4.0937E - 03 2.432E—02 6.057E—03 k2o aThis is the difference between the observed and calculated values multiplied by the weighting factor.

Interspecies and Dose-Route Extrapolations Curtis C. Travis INTRODUCTION Risk assessment is a procedure used to synthesize all available data and the best scientific judgment to estimate the risk associated with human exposure to chemicals. Because of gaps in our current scientific under- standing of the cancer-causing process, risk assessment requires the use of a series of judgmental decisions on unresolved issues. The major as- sumptions arise from the necessity to extrapolate experimental results (1) across species from rats or mice to humans and (2) from the high- dose regions to which animals are exposed in the laboratory to the low- dose regions to which humans are exposed in the environment, and (3) across routes of administration. There is a growing awareness of the need for an evaluation of the scientific bases for the assumptions used in the risk assessment process. Pharmacokinetics provides a tool for such an eval- uation. Pharmacokinetics is the study of the absorption, distribution, metabo- lism, and elimination of chemicals in man and animals. An effective approach for interpreting empirical data relating to pharrnacokinetics is the development of predictive physiologically based pharmacokinetic mod- els. These models utilize actual physiological parameters of the experi- The submitted manuscript has been authored by a contractor of the U.S. Government under contract DE-AC05-840R21400. Accordingly, the U.S. Government retains a nonexclusive roy- alty-free license to publish or reproduce the published form of this contribution. or allow others to do so, for U.S. Government purposes. 208

Interspecies and Dose-Route Extrapolations 209 mental animals such as breathing rates, blood flow rates, tissues volumes, etc., to describe the pharmacokinetic process. These physiological param- eters are coupled with chemical-specific parameters such as blood/gas coefficients, tissue/blood partition coefficients, and metabolic constants to predict the dynamics of a compound's movement through an animal system. The models have the ability to relate exposure concentrations quantitatively to organ concentrations over a range of exposure intervals. As such, the models allow for prediction of the relationship between inhaled concentrations of a chemical and the concentration found in target tissues. A chief advantage of the physiologically based model is that by simply changing the physiological parameters, the same model can be utilized to describe the dynamics of chemical transport and metabolism in mice, rats, and humans. DESCRIPTION OF THE MODEL The pharmacokinetic model used in the present study (Figure 1 and Table 1) is patterned after that developed by Dedrick (1973) and modified by Ramsey and Andersen (1984) for prediction of the behavior of styrene inhalation exposure in humans from behavior observed in rats. The chem- icals used as an illustration in this particular study are tetrachloroethylene (perchloroethylene, Perc, PCE), which is used extensively as a dry-clean- ing solvent and methylene chloride (dichloromethane, DCM), which is a solvent with wide use in the food industry. Details of the model devel- opment are given in Word et al. (in press). The tissue groups include (1) vessel-rich organs such as brain, kidney, and viscera; (2) vessel-poor organs such as the muscle and skin; (3) slowly perfused fat tissue; and (4) organs with a high capacity to metabolize (principally liver). The model is described mathematically by a set of differential equations which calculate the rate of change of the amount of chemical in each compartment. Metabolism of PCE and DCM, which occurs chiefly in the liver, is described by a combination of a linear metabolic component and a Michaelis-Menten component describing sat- urable metabolism. In the modeling of the uptake, distribution, and elimination of PCE and DCM, we used the 0.75 exponential power of the body weight (BOO) to scale cardiac output and ventilation rate within species. The scaling for- mulas are as follows: Call = Fail` BW075, and Qb = Qb~ BW0 7s, (1) (2) where the allometric constants Qa~vc and Qb~ are found in Table 2. The blood flow to a given tissue group is obtained from the total blood flow

21~ 0 CURTIS C. TRAVIS M ETABO L I TES ( L I N EA R PATHWAY ) Qalv ~ ALVEOLAR _ C i n h SPAC E Qb LUNG C BLOOD ven FAT T I SSU E Cvf G ROUP V ESSE L POO R Cv p G R OU P VESSE L R ICH Cv r G ROU P LIVE R METABO LIZ I NG Cvl TISSUE G ROUP l l ~J I Vmax, Km M ETABO LITES (M ICHAE LIS-MENTEN ) Qalv Calv Qb Cart Qf Cart Q. I Cart Ql Cart FIGURE 1 Diagram of the pharmacokinetic model used to simulate the behavior of inhaled PCE. The symbols are defined in Table 1, and the parameters used to describe the model are described in Table 3.

Interspecies and Dose-Route Extrapolations 211 TABLE 1 Abbreviations and Symbols Used in Describing a Physiologically Based Pharmacokinetic Model Am Qi Vj Cj phi Do Alveolar ventilation rate (liters air/in) Concentration in inhaled air (mg/liter air) Concentration in alveolar air (mg/liter air) Blood/air partition coefficient (liters air/liters blood) Cardiac output (liters blood/in) Concentration in arsenal blood (mg/liter blood) Concentration in mixed venous blood (mg/liter blood) Michaelis-Menten metabolism rate (mg/h) Michaelis constant (mg/liter blood) Linear metabolism rate (h ~ I) Amount metabolized in the liver (mg) Subscripts (i) for tissue groups or compartments I Liver (metabolizing tissue group) f Fat tissue group r Vessel-rich tissue group p Vessel-poor tissue group Blood flow rate to tissue group (liters blood/in) Volume of tissue group i (liters) Concentration in tissue group i (mg/liter) Amount in tissue group i (mg) Concentration in venous blood leaving tissue group i (mg/liter blood) Tissue/blood partition coefficient for tissue i (liters blood/liter)) Gavage or oral rate constant (h ~ ') Total quantity of PCE absorbed via gavage route (mg) Qb by multiplying by the fraction of blood flowing to that tissue. The fractions are given in Table 2. Likewise, the volume of a tissue group is determined by multiplying the body weight by the volume fractions given in Table 2. The bone and cartilage volume (9%) is ignored in the model. The tissue/blood and tissue/air partition coefficients used in our model were made available to us by M. E. Andersen (from D. M. Hetrick, Oak Ridge National Laboratory, Oak Ridge, Tenn., personal communication). The partition coefficients were measured by using a vial equilibration technique (Sato and Nakajima, 1979) in which the chemical was added to a closed vial containing blood or tissue, and the partitioning was de- termined by estimating the amount that disappeared from the head space after equilibration at 37°C (Gargas et al., 1986~. The tissue/blood partition coefficients used in the model (Table 3) are obtained by dividing the tissue/air partition coefficients by the blood/air partition coefficients. Gar- gas et al. (1986) have also proposed that the metabolic parameters Vma,` and Kf can be scaled using body weight. The Michaelis-Menten metabolic

212 CURTIS C. TRAVIS TABLE 2 Coefficients for Scaling Formulas Rats Humans Allometr~c constants Ventilation (Qa~vc; liters/in) 14 14 Cardiac output (Qbc; liters/in) 14 14 Tissue volume fractions Liver 0.04 0.04 Fat 0.07 0.19 vRGa 0.05 °.°S VPGb 0.75 0.63 Blood flow fractions Liver 0.25 0.25 Fat 0.O9 0.O9 vRGa 0.5 1 0.5 1 vpGb 0.15 0.15 aVessel-nch group. bVessel-poor group. rate Vma,; is assumed to be proportional to surface area (Andersen et al. 1984): Vma,` = Vm~,~c BW0 7. The linear metabolic rate constant is scaled as follows: Kf = KfC BW -0 3. (3) (4) Values for the coefficients Vma,~c and KfC and for the Michaelis-Menten constant Km are shown in Table 4. INTERSPECIES EXTRAPOLATION In this section, it will be determined whether the metabolic scaling factors described above allow for both interspecies extrapolation and dose- route extrapolation. To accomplish this, we will first determine metabolic parameters so that model predictions will reproduce rat data published by Pegg et al. (1979). Then, using scaling parameters, we will attempt to reproduce rat ingestion data published by Pegg et al. (1979) and human data published by Fernandez et al. (1976). Rat Inhalation Adult male Sprague-Dawley rats weighing 250 g were exposed to ~4C- labeled PCE by inhalation for a duration of 6 h in experiments conducted

Interspecies and Dose-Route Extrapolations 213 TABLE 3 Physiological and Biochemical Parameters Used in Describing the Behavior of PCE in the Pharmacokinetic Model of Figure 1 Parameter Rat Human Body weight (kg) BW 0.25 JO oa Alveolar ventilation (liters air/in) Qalv 5.02 325.0 Blood flow rates (liters blood/in) Total blood flow rate Qb 5.02 325.0 Blood flow rate in liver I 1.26 81.3 Blood flow rate in fat Q5 0.45 29.3 Blood flow rate in vessel-nch tissues Or 2.56 165.7 Blood flow rate in vessel-poor tissues Qp 0.75 48.7 Tissue group volumes (liters) Volume in liver Vat 0.0100 2.8 Volume in fat Vf 0.0175 13.3 Volume in vessel-nch tissues Vr 0.0125 3.5 Volume in vessel-poor tissues Vp 0.1875 44.1 Blood/air partition coefficient Ab 18.9 10.3 Tissue/air partition coefficients Liver/air partition coefficient Airs 70 3 70 3 Fat/air partition coefficient Af/a 2,060.0 2,060.0 Vessel-rich/air partition coefficient Aria 70 3 70 3 Vessel-poor/air partition coefficient Apia 20.0 20.0 aWith the exception of Fernandez et al. (1976) (83 kg). by Pegg et al. (19791. In the 72 h following exposure to 10 ppm, metab- olism accounted for 20% of the total radioactivity recovered, while un- changed PCE in expired air accounted for 70%. Pulmonary elimination of PCE had a half-life of about 7 h. The biological parameters and partition coefficients used to model the empirical data of Pegg et al. (1979) are presented in Table 3. We used Andersen's values for the Atla (tissue/air) partition coefficients, except for fat. It was necessary to increase the fat/air partition coefficient from Andersen's measured value of 1,638 to 2,060 to account for increased alveolar concentrations of PCE after exposure. Since empirical values for the metabolic parameters are not available, we determined values for these parameters which produced the best fits with the inhalation data from Pegg et al. (19791. Metabolic scaling coefficients which produced the best fit with the empirical data were Vma,`— 0.068 mg/h, Km = 0.3 mg/liter blood, and Kf = 2.73/h. In Figure 2, model predictions are compared with the empirical data of Pegg et al. (19791. Figure 2 shows the percentage of PCE recovered in expired air in rats (for a number of different time intervals) following an exposure to 10 ppm of PCE for 6 h. The vertical bars represent the range

214 CURTIS C. TRAVIS 20 — J ~ ~ 15 ~ I_ O Z cat tl] ~ _ flu ~ ~ tr 1 0 O O LU c, ~ He LU con ~ 5 x UJ UJ o 1 10 ppm PCE EXPOSURE IN RATS FOR 6 h - · MODEL PREDICTIONS T EXPERIMENTAL DATA l WITH ERROR BAR I I l · ~ 0 8 16 24 32 40 48 56 64 72 POSTEXPOSURE TIME (fur) FIGURE 2 Percentage of PCE expired by rats following expsure to to ppm in air for 6 h. The percentage plotted at 8 h represented the percentage expired from 6 to 8 h, etc. Note that the time intervals are longer after 8 h. Symbols: ., model predictions; bars, range of experimental data.Source: Pegg et al. (1979). Ot the empirical data. The predictions and empirical data are integrated from the beginning of each time interval to the time at which the points and bars are drawn. Note that the length of the time intervals varies, causing the data points to be nonmonotonic. The model predicted that 68% of the body burden of PCE would be recovered in expired air during the 72 h after exposure, which is in good agreement with the experimen- tally determined value of 70%. Rat Ingestion Adult male Sprague-Dawley rats weighting 250 g were orally admin- istered ~4C-labeled PCE in corn oil in experiments conducted by Pegg et

Interspecies and Dose-Route Extrapolations 215 al. (19791. In the 72 h following exposure to 500 mg/kg, metabolism accounted for 5% of total radioactivity recovered, while unchanged PCE in expired air accounted for 90%. There was no significant difference in the elimination half-life (approximately 7 h) with dose or route of ad- ministration. In Figure 3, model predictions of expired air concentrations are com- pared with the empirical data of Pegg et al. (19791. The model, based on parameters obtained from the inhalation study, slightly overpredicted the rate of elimination of PCE. This can be attributed to the effects of corn oil as a carrier vehicle on the pharmacokinetics of PCE. Withey et al. (1983), Withey (1984), and Angelo et al. (1986) have shown that an oil carrier results in a slower elimination pattern for dichloromethane (meth- ylene chloride). 10 J it_ Icy 6 1 0.] tY cat x ~ 0.01 cat 0.001 ^L r~,1 L1 \t] ~ \ lo- \ A ~ \ ~ . ~ \ \ ~1 l ll [a -1 ~ 1 ~~ 0 720 1440 2160 2880 3600 4320 TIME (/~IN) PERC RATS GAVAGE (139~1N) FIGURE 3 Concentration of PCE expired by rats following oral administration of 500 mg of i4C-labeled PCE per kg in corn oil for 6 h. Symbols: solid line, model predictions; O. range of experimental data.

2 ~ 6 CU RTIS C. TRAVIS Human Inhalation Fernandez et al. (1976) exposed 24 subjects to concentrations of 100 ppm of PCE for 1 to 8 h. During the first hours after exposure, the concentration of PCE decreased rapidly in alveolar air; however, more than 2 weeks was necessary to eliminate the PCE retained following an exposure of 100 ppm for 8 h. The biological parameters and partition coefficients used to model the data of Fernandez et al. (1976) are also presented in Table 3. A body weight of 70 kg was used in all model calculations except for those of Fernandez et al., in which body weight was 83 kg. Metabolic parameters for humans were obtained using the scaling for- mulas and the Vm~,,~c and Kfc determined for rats (Table 4~. Figure 4 is the graph of the data of Fernandez et al. (1976), showing alveolar concentrations resulting from a 100-ppm exposure to PCE for various exposure durations, versus our model results. The model results were calculated assuming 19% fat, and the agreement between the model and the data is very good. DOSE-ROUTE EXTRAPOLATION The route of exposure to organic chemicals can significantly affect the quantity of a chemical that reaches a particular target tissue (Angelo et al., 1986; National Research Council, 19861. Pharmacokinetic models provide a tool to quantitatively evaluate the effect of route of administration on dose to the target tissue. This effect must be evaluated on a chemical by chemical basis. The standard procedures for calculating applied dose following inhal- ation or ingestion exposures are as follows. For inhalation, the applied dose (in milligrams) is the product of the air concentration (in milligrams per liter) of the chemical, the animal breathing rate (in liters per minute), and the duration of exposure (in minutes). For drinking water ingestion, the applied dose (in milligrams) is the product of the water concentration (in milligrams per liter) of the chemical, the animal drinking water rate TABLE 4 Metabolic Parameters Used in Describing the Behavior of PCE Parameter Rat Human Body weight (kg) BW 0.25 70.0 Michaelis-Menten metabolic rate (mg/h) imp 0.068 3.5 Michaelis constant (mg/liter blood) Km 0~3 0 3 Linear metabolic rate (h- ') Ks 2.73 0.0

Interspecies and Dose-Route Extrapolations 217 70 60 50 40 Q - ~0 30 20 10 o 1 1 1 1 1 1 1 1 l hi I Xe - ,2 hr \4 hr · ·W \ ALVEOLAR CONCENTRATIONS (ppm) FOR 100 ppm PCE · EXPOSURE IN HUMANS — MODEL PREDICTIONS ·l : \ · ·` —.\8 hr `: ·\ · - · EXPER I M ENTAL DATA O F FER NAN D EZ et a/. (1976) - - - 0 2 4 6 8 10 12 14 16 TIME (h) FIGURE 4 PCE concentrations in alveolar air of humans after an exposure concentration of 100 ppm for periods of 1-, 2-, 4-, and 8-h durations. Symbols: solid line, model predictions; ., experimental data. (in liters per minute), and the duration of exposure (in minutes). In both of these formulas, 100% absorption into the body has been assumed, a standard assumption in risk analysis when data to the contrary are lacking. The pharmacokinetic model developed by Andersen et al. (1987) was applied to investigate the dose to target tissues following both inhalation and oral administration of DCM in mice. Figure 5 shows model predictions of total metabolized dose to the mouse liver over a 24-h period. The inhalation exposure was for a duration of 6 h and the ingestion exposure was for 24 h. The applied dose for inhalation exposures was calculated using the alveolar ventilation rate as opposed to the minute volume. Figure 5 indicates fairly good agreement in the effective liver dose for the two routes of administration. The largest difference is a factor of three and occurs in the 500- to 10,000-mg/kg applied dose range. Figure 6 shows model predictions of total metabolized dose to the mouse lung following inhalation and oral administration. As can be seen, the effective lung dose shows greater dependence on the route of administra- tion. The largest difference is a factor of 5 and, in contrast to the liver, occurs in the low applied dose range.

218 CURTIS C. TRAVIS 1000 100 10 - L~ En o C) hJ hJ LO 0.1 0.01 0.001 DRINKING WATER INGESTION / // LIVER // // // // / INH>ATION 1 1 1111111 1 1 1111111 1 1 1111111 1 1 1111' // .....n lll l l l l l l 10 100 1000 10000 100000 1000000 APPLIED DOSE (MG/KG) FIGURE 5 Effective dose to the mouse liver from inhalation and oral ingestion of DCM. DISCUSSION It has been shown that it is possible to extrapolate pharmacokinetic responses to the chemical PCE by scaling of physiological parameters. Simple scaling of metabolic parameters, originally adjusted to account for the effects in rats, also accounted for the responses in human subjects to doses of PCE administered by inhalation. It has also been shown that it is possible to use pharmacokinetic models to investigate the effect of route of administration on dose to target tissue. Since pharmacokinetic models allow for a quantitative extrapolation of exposure data across species and between routes of administration, they provide a tool to quantitatively evaluate assumptions currently used in the risk assessment process. While the present paper dose not specifically evaluate current assumptions, it does demonstrate that such an evaluation is presently possible. 1

Interspecies and Dose-Route Extrapolations 219 1000 100 10 u, o 1 0.1 0.01 0.001 LUNG DRINKING WATER INGESTION ,.~ // IN~ATION // 1 1 1111111 1 1 1111111 1 1 1111111 1 1 1111111 1 1 111111 1 10 100 1000 10000 100000 1000000 APPLIED DOSE (MG/KG) FIGURE 6 Effective dose to the mouse lung from inhalation and oral ingestion of DCM. REFERENCES Andersen, M. E., R. L. Archer, H. J. Clewell III, and M. G. MacNaughton. 1987. A physiological model in the rat. The Toxicologist 4:111. Andersen, M. E., H. J. Clewell III, M. L. Gargas, F. A. Smith, and R. H. Reitz. 1987. Physiologically based pharmacokinetics and the risk assessment process for methylene chloride. Toxicol. Appl. Pharmacol. 87:185. Angelo, M. J., A. B. Pritchard, D. R. Hawkins, A. R. Waller, and A. Roberts. 1986. The pharmacokinetics of dichloromethane I: Disposition in B6C3F~ mice following in- travenous and oral administrations. Food Chem. Toxicol. 24:965. Dedrick, R. L. 1973. Animal scale-up. J. Pharmacol. Biopharm. 1:435. Fernandez, J., G. Guberan, and J. Caperos. 1976. Experimental human exposures to tetrachloroethylene vapor and elimination in breath after inhalation. Am. Indust. Hyg. Assoc. J. 37:143. Gargas, M. L., H. J. Clewell III, and M. E. Andersen. 1986. Metabolism of inhaled dihalomethanes in vivo: Differentiation of kinetic constants for two independent pathways. Toxicol. Appl. Pharmacol. 82:211. National Research Council. 1986. Drinking Water and Health, Vol. 6. Washington, D.C. National Academy Press.

220 CURTIS C. TRAVIS Pegg, D. G., J. A. Zempel, W. B. Braun, and P. G. Watanabe. 1979. Disposition of tetrachloro (TIC) ethylene following oral and inhalation exposure in rats. Toxicol. Appl. Pharmacol. 51:465. Ramsey, J. C., and M. E. Andersen. 1984. A physiologically based description of the inhalation pharmacokinetics of styrene in rats and humans. Toxicol. Appl. Pharmacol. 73:159. Sato, A., and T. Nakajima. 1979. Partition coefficients of some aromatic hydrocarbons and ketones in water, blood and oil. Br. J. Indust. Med. 36:231. Ward, R. C., C. C. Travis, D. M. Hetrick, M. E. Andersen, and M. L. Gargas. In press. Pharmacokinetics of tetrachloroethylene. Toxicol. Appl. Pharmacol. Withey, J. R. 1984. Classical pharmacokinetics of methylene chloride-oral administration. In Proceedings of the Food Solvents Workshop I: Methylene Chloride. Washington, D.C.: Nutrition Foundation. Withey, J. R., B. T. Collins, and P. G. Collins. 1983. Effect of vehicle on the pharma- cokinetics and uptake of four halogenated hydrocarbons from the gastrointestinal tract of the rat. J. Appl. Toxicol. 3:249.

Carcinogen-DNA Adduces as a Measure of Biological Dose for Risk Analysis of Carcinogenic Data Marshall W. Anderson One of the major problems confronting the regulatory agencies today is the extrapolation of high-dose toxicology data in animals for the as- signment of potential risk in the human population. Both low-dose and species-to-species extrapolation of toxicology data are needed for proper risk assessment of human exposure to chemicals. It is recognized that dose-response relationships for the toxicological response under considerations are of great value in constructing extrap- olation procedures. Although it is not practical in most cases to measure toxicological responses, neoplasia, immunosuppression, etc., at relevant environmental or occupational exposure levels, a measurement of dose other than exposure dose may enhance the sensitivity of the risk analysis of high-dose toxicological data. The choice of the biological dose should be based on the mechanisms involved in the toxicological response under consideration. This report will discuss the utility of using carcinogen- DNA adduct levels as a measure of the biological dose in the risk analysis of carcinogenic data. EVIDENCE FOR USE OF DNA ADDUCTS AS A MEASURE OF BIOLOGICAL DOSE Considerable evidence has indicated that many mutagens and carcin- ogens react with cellular DNA either directly or following metabolic for- mation of reactive products. If DNA replication proceeds on such a modified template before altered bases or nucleotides are removed by enzymatic 221

222 MARSHALL W. ANDERSON repair processes, the mutations can be genetically fixed. Thus, the extent of promutagenic damage induced by environmental chemicals and the capacity of the cell to repair such damage may be important factors in both initiation of malignant transformation and tissue specificity of many carcinogens. Various studies, both in vivo and in vitro, of carcinogen-DNA adducts in a known target tissue are a good measure of a biological dose for initiation of neoplasia. 1. Data on the mutagenicity of many chemical carcinogens in vitro implies that mutation results from the cells' attempt to deal with unexcised DNA adducts at the time of replication (Fatal et al., 1981; McCormick and Maher, 1985; Newbold et al., 1979; van Zeeland et al., 1985; Yang et al., 19801. A similar conclusion was reached by Russell et al. (1982a,b) in a study of ethylnitrosourea-induced transmitted mutations in mice. 2. For equal levels of O6-ethylguanine in the DNA of V798 cells and in testicular DNA from male mice treated with ethylnitrosourea, the fre- quency of mutation induction in V79 cells in vitro and in the specific- locus assay in the mouse were very similar (van Zeeland et al., 19851. 3. A positive correlation between the carcinogenicity of a series of polycyclic aromatic hydrocarbons and their extent of reaction with DNA has been observed (Brookes and Lawley, 1964; Goshman and Hiedel- berger, 1967; Huberman and Sachs, 19671. The binding of 3-propiolactone and other alkylating agents to DNA correlates with their tumor-initiating potency (Colburn and Boutwell, 19681. 4. Several investigators have examined the effect of inhibitors of car- cinogenesis on the formation of adducts between carcinogen metabolites and DNA. In general, the degree of inhibition of tumor induction correlates with the degree of inhibition of adduct formation. For example, benzoVa~pyrene diol epoxide (BPDE)-DNA adduct formation and benzoka~pyrene-induced neoplasia in mouse lung, forestomach, or skin were inhibited to the same degree by butylated hydroxyanisole (Anderson et al., 1981, 1985; Adriaenssens et al., 1983; Ioannou et al., 1982), inducers of aryl hydrocarbon hydroxylase (Anderson and Bend, 1983; Cohen et al., 1979; Ioannou et al., 1982; Wilson et al., 1981), and a- angelicalactone (Ioannou et al., 19821. Treatment with 2,3,7,8-tetrachlor- odibenzo-p-dioxin (TCDD) completely inhibited papilloma formation in- duced by 7,12-dimethylbenzEal-anthracine (DMBA) and DMBA adduct formation (Cohen et al., 1 9791. 5. Swann et al. (1980) showed that changes in the incidence of di- methylnitrosamine (DMN)-induced kidney tumors produced by changes in the diet and by treatment with benzofaypyrene correspond to the changes

Carcinogen-DNA Adducts in Risk Analysis 223 these treatments produce in the alkylation of the target tissue DNA by DMN. 6. Janss and Benn (1978) found a correlation between the amount of DMBA bound to DNA and the incidence of mammary tumors in rats of different ages. 7. Results from Richardson et al. (1986) demonstrated that interlobe differences in the incidence of diethylnitrosamine-induced hepatocellular carcinomas are due in part to differences in DNA alkylation and cell replication. Thus, the heterogenous tumor response within an organ is quantitatively related to the molecular dosimetry of DNA adducts. 8. Several studies have shown that DNA adducts preferentially accu- mulate in target cell DNA during continuous exposure to carcinogen. When rats are continuously exposed to hepatocarcinogenic regimens of the al- kylating agents dimethyluitrosamine or dimethylhydrazine, 06-methyl- guanine accumulates in the target nonparenchymal cells but not in hepatocytes (Bedell et al., 1982; Lindamood et al., 19821. The Clara cell, although accounting for only 1% pulmonary cells in the lungs of rats, was found to possess a much higher concentration of 06-methylguanine than the other lung cell types, particularly at lower doses (Belinsky et al., 19871. The Clara cell is the putative progenitor cell for lung tumors induced by 4-(N-methyl-N-nitros amino)-1-~3-pyridyl)-1-butanone (NNK). 9. Recent studies on oncogene activation in rodent tumors have shown that the activating mutations observed in ras genes are consistent with the known DNA adduct patterns of the carcinogens (Barbacid, in press). In general, these observations demonstrate a correlation between spe- cific carcinogen-DNA adduct levels in the target tissue and tumor response. Tissue concentration of carcinogen-DNA adducts, however, may not ex- plain the differences in organ susceptibility and species susceptibility to chemically induced neoplasia. Other aspects in the multistep development of neoplasia may be required to explain these differences. In any case, this does not detract from the use of adduct levels in the target tissue as a measure of biological dose of the carcinogen. The abovementioned results strongly imply that it is biologically more meaningful to relate tumor response to concentrations of specific DNA adducts in the target tissue than it is to relate tumor response to the administered dose of the chemical (Hoer et al., 19831. FACTORS TO CONSIDER IN CONSTRUCTING A MEASURE OF BIOLOGICAL DOSE FROM CARCINOGEN-DNA ADDUCT LEVELS Carcinogen-DNA adducts only represent potential promutagenic le- sions. Mutations are genetically fixed only if DNA replication proceeds

224 MARSHALL W. ANDERSON on the modified template before altered bases or nucleotides are removed by DNA repair processes. Thus, the extent of cell replication must be considered in addition to the accumulation of specific DNA adducts. Swenberg et al. (1983, 1985) suggest that the product of (cell replication) x (the concentration of carcinogen-DNA adduct) x (the number of cells at risk) is one possible measure of biological dose. They used this cal- culation of initiation index to explain the cell specificity in hepatocarcin- ogenesis during continuous exposure of rats to 1,2-dimethylhydrazine (Bedell et al., 1982) and mice to dimethylnitrosamine (Lindamood et al., 19821. DETERMINATION OF CARCINOGEN-DNA ADDUCT ACCUMULATION The accumulation of carcinogen-DNA adducts in a cell is the difference between adduct formation and removal of adducts by DNA repair pro- cesses. The dose- and time-dependent accumulation of adducts is required to construct a measure of biological dose. Adduct levels can be determined from direct measurements or, at least in theory, calculated a priori from physiologically based pharmacokinetic models. As shown by several pa- pers in this volume, this type of model can be used to predict blood or tissue levels of the parent compound or metabolites under a variety of experimental conditions. The complexity of adduct accumulation, how- ever, suggests that this same approach may not be feasible to previously predicted adduct levels. Adduct levels comprise a very small percentage of the total administered dose, as small as 10-5% of the total dose (Ad- riaenssens et al., 19831. Moreover, it may be necessary to determine levels of adducts in selective cell types to obtain adequate measures of biological dose, as illustrated by several reports (Bedell et al., 1982; Belinsky et al., 1987; Lindamood et al., 19821. Prediction of small quantities in individual cell types would probably require very accurate measurements of a large number of metabolic rate constants and, in addition, determination of enzymatic parameters that could determine repair rates of specific carcin- ogen-DNA adducts. At present, it is probably more feasible to measure directly DNA adduct levels. Sensitive and accurate procedures have been developed to measure adducts at levels of 1 modified base per 109 normal nucleotides and potentially at even lower levels of DNA damage (Adamkiewicz et al., 1985; Baan et al., 1985; Balhorn et al., 1985; Fisher et al., 1985; Randerath et al., 1985~. Even though the types of pharma- cokinetic models constructed to predict blood and tissue levels of parent compounds or metabolites may not predict accumulation of adduct in selective cell types, it may be possible to construct similar types of kinetic models for adduct accumulation.

Carcinogen-DNA Adducts in Risk Analysis 225 Chemicals may damage DNA by mechanisms other than forming stable DNA adducts. For example, chemically induced free radical formation may cause DNA strand breaks and/or DNA-DNA or DNA-protein cross- links. Chemically, unstable adducts may generate apurinic/apyrimidinic sites which can be mutagenic (Vousden et al., 1986~. More studies are required to evaluate the potential risk from these types of DNA lesions as well as to develop sensitive assays to measure these lesions. DOSE-RESPONSE RELATIONSHIPS - Dose-response relationships for carcinogen-DNA adducts have been determined for several chemicals. A plot of adduct levels divided by dose versus dose is one way to represent dose-response relationships, and this representation is especially useful for consideration of low-dose extrap- olation. The three types of curves illustrated in Figure 1 have been ob- sel~ved. For relationships like that depicted in curve A of Figure 1 (Ioannou et al., 1982), the exposure dose and adduct levels are interchangeable \ _~ B/ DOSE FIGURE 1 Relationships between exposure dose and ad- duct levels. See the text for an explanation of the three curves. I'

226 MARSHALL W. ANDERSON because a linear relationship exists. For relationships like that depicted in curve B of Figure 1 (Hoer et al., 1983; Ioannou et al., 1982), the adduct levels as a percentage of dose decreases as the dose decreases, and thus, exposure dose would overestimate low-dose risk. For relationships like that depicted in curve C of Figure 1 the adduct levels as a percentage of dose increases as the dose increases (Belinsky et al., 1987), and thus, exposure dose would underestimate low-dose risk. For continuous exposure to carcinogen, adduct levels must be measured at various time points for each dose. The biological dose would then be related to the area under the time curve for adduct levels. Again, dose- dependent cell replication must be considered in constructing the biological dose. Other modes of action of the chemical related to tumor induction should also be considered when information is available. In any case, for purposes of low-dose extrapolation, adduct levels as a measure of bio- logical dose is superior to exposure levels since for situations like curve C, low-dose risk will not be underestimated. REFERENCES Adamkiewicz, J., G. Eberle, N. Huh, P. Nehls, and M. Rajewsky. 1985. Quantitation and visualization of alkyl deoxynucleosides in the DNA of mammalian cells by mono- clonal antibodies. Environ. Health Perspect. 62:49-55. Adriaenssens, P. I., C. M. White, and M. W. Anderson. 1983. Dose-response relationships for the binding of benzo(a)pyrene metabolites to DNA and protein in lung, liver, and forestomach of control and butylated hydroxyanisole-treated mice. Cancer Res. 43:3712- 3719. Anderson, M. W., and J. R. Bend. 1983. In viva metabolism of benzo-(a)-pyrene: For- mation and disappearance of BB-metabolite-DNA adducts and extrahepatic tissues vs. liver. Pp. 459-467 in J. Ragdstrom, J. Montelius, and M. Bengtsson, eds. Extrahepatic Drug Metabolism and Chemical Carcinogenesis. New York: Elsevier Science Publishers. Anderson, M. W., M. Boroujerdi, and A. G. E. Wilson. 1981. Inhibition in viva of the formation of adducts between metabolites of benzo(a)pyrene and DNA by butylated hydroxyanisole. Cancer Res. 41 :4309-4315. Anderson, M. W., P. I. Adriaenssens, C. M. White, Y. M. Ioannou, and A. G. E. Wilson. 1985. Effect of the antioxidant butylated hydroxyanisole on in viva formation of benzo- (a)-pyrene metabolite-ONA adducts. P. 241 in J. W. Finley and D. E. Schwass, eds. Xenobiotic Metabolism: Nutritional Effects. ACS Symposium Series No. 277. Wash- ington, D.C.: American Chemical Society. Baan, R. A., O. B. Zaalberg, and A. M. J. Fichtinger-Schepman, M. A. Muysken-Schoen, M. J. Lansbergen, and P. H. M. Lohman. 1985. Use of monoclonal and polyclonal antibodies against DNA adducts for the detection of DNA lesions in isolated DNA in single cells. Environ. Health Perspect. 62:81-88. Balhorn, R., J. A. Mazrimas, and M. Corzett. 1985. Application of HPLC to the isolation of molecular targets in dosimetry studies. Environ. Health Perspect. 62:73-79. Barbacid, M. In press. Annul Rev. Biochem. Bedell, M. A., J. G. Lewis, K. C. Billings, and J. A. Swenberg. 1982. Cell specificity in hepatocarcinogenesis: Preferential accumulation of 06-methylguanine in target cell

Carcinogen-DNA Adducts in Risk Analysis 227 DNA during continuous exposure to rats to 1,2-dimethylhydrazine. Cancer Res. 42:3079- 3083. Belinsky, S. A., C. M. White, J. A. Boucheron, F. C. Richardson, J. A. Swenberg, and M. W. Anderson. 1987. Cell selective alkylation of DNA in rat lung following low dose exposure to the tobacco specific carcinogen 4-(N-methyl-N-nitrosamino)-1-(3-pyridyl)- 1-butanone. Cancer Res. 47: 1143- 1148. Brookes, P., and P. D. Lawley. 1964. Evidence for the binding of polynuclear aromatic hydrocarbons to nucleic acid of mouse skin; relation between carcinogenic power of hydrocarbons and their binding to deoxyribonucleic acid. Nature 202:781-784. Cohen, G. M., W. M. Bracken, R. P. Iyer, D. L. Berry, J. K. Selkirk, and T. J. Slaga. 1979. Anticarcinogenic effects of 2,3,7,8-tetrachlorodibenzo-p-dioxin on benzo(a)pyrene and 7,12-dimethylbenz(a)anthracene tumor initiation and its relationship to DNA binding. ~ Cancer Res. 39:4027-4033. Colburn, N. H., and R. K. Boutwell. 1968. The binding of beta-propiolactone and some related alkylating agents to DNA, RNA, and protein of mouse skin; relation between tumor-initiating power of alkylating agents and their binding to DNA. Cancer Res. 28:653-660. Fahl, W. E., D. G. Scarpelli, and K. Gill. 1981. Relationship between benzo(a)pyrene- induced DNA base modification and frequency of reverse mutations in mutant strains of Salmonella typhimurium. Cancer Res. 41 :3400-3406. Fisher, D. H., J. Adams, and R. W. Giese. 1985. Trace derivatization of cytosine with pentafluorobenzoyl chloride and dimethyl sulfate. Environ. Health Perspect. 62:67-71. Goshman, L. M., and C. Heidelberger. 1967. Binding of tritium-labeled polycyclic hy- drocarbons to DNA of mouse skin. Cancer Res. 27:1678-1688. Hoel, D. G., N. L. Kaplan, and M. W. Anderson. 1983. Implication of nonlinear kinetics on risk estimation in carcinogenesis. Science 219:1032-1037. Huberman, E., and L. Sachs. 1967. DNA binding and its relationship to carcinogenesis by different polycyclic hydrocarbons. Int. J. Cancer 19: 122- 127. Ioannou, Y. M., A. G. E. Wilson, and M. W. Anderson. 1982. Effect of butylated hy- droxyanisole, alpha-angelica lactone, and beta-naphthoflavone on benzo(alpha)pyrene: DNA adduct formation in vivo in the forestomach, lung, and liver of mice. Cancer Res. 42:1199-1204. Janss, D. H., and T. L. Benn. 1978. Age-related modification of 7,12-demethyl- benz(a)anthracine binding to rat mammary gland DNA. J. Natl. Cancer. Inst. 60:173- 177. Lindamood, C., M. A. Bedell, K. C. Billings, and J. A. Swenberg. 1982. Alkylation and de novo synthesis of liver cell DNA from C3H mice during continuous dimethylnitro- samine exposure. Cancer Res. 42:4153-4157. McCormick, J. J., and V. M. Maher. 1985. Cytotoxic and mutagenic effects of specific carcinogen-DNA adducts in diploid human fibroblasts. Environ. Health Perspect. 62:145- 155. Newbold, R. F., P. Brookes, and R. G. Harvey. 1979. A quantitative comparison of the mutagenicity of carcinogenic polycyclic hydrocarbon derivatives in cultured mammalian cells. Int. J. Cancer 24:203-209. Randerath, K., E. Randerath, H. P. Agrawal, R. C. Gupta, M. E Schurdak, and M. V. Reddy. 1985. Postlabeling methods for carcinogen-DNA adduct analysis. Environ. Health Perspect. 62:57-65. Richardson, F. C., J. A. Boucheron, M. C. Dyroff, J. A. Popp, and J. A. Swenberg. 1986. Biochemical and morphologic studies of heterogeneous lobe responses in hepa- tocarcinogenesis. Carcinogenesis 7:247-251.

228 MARSHALL W. ANDERSON Russell, W. L., P. R. Hunsicker, D. A. Carpenter, C. V. Cornett, and G. M. Guinn. 1982a. Effect of dose fractionation on the ethylnitrosourea induction of specific-locus mutations in mouse spermatogonia. Proc. Natl. Acad. Sci. USA 79:3592-3593. Russell, W. L., P. R. Hunsicker, G. D. Raymer, M. H. Steele, K. F. Stelzner, and H. M. Thompson. 1982b. Dose-response curve for ethylnitrosourea-induced specific-locus mu- tations in mouse spermatogonia. Proc. Natl. Acad. Sci. USA 79:3589-3591. Swann, P. F., D. G. Kaufman, P. N. Magee, and R. Mace. 1980. Induction of kidney tumors by a single dose of dimethylnitrosamine: Dose response and influence of diet and benzo(a)pyrene pretreatment. Br. J. Cancer 41:285-294. Swenberg, J. A., D. E. Rickert, B. L. Baranyi, and J. I. Goodman. 1983. Cell specificity in DNA binding and repair of chemical carcinogens. Environ. Health Perspect. 49:155- 163. Swenberg, J. A., F. C. Richardson, J. A. Boucheron, and M. C. Dyroff. 1985. Relation- ships between DNA adduct formation and carinogenesis. Environ. Health Perspect. 62: 177- 183. van Zeeland, A. A., G. R. Mohn, A. Neuhauser-Klaus, and V. H. Ehling. 1985. Quan- titative comparison of genetic effects of ethylating agents on the basis of DNA adduct formation. Use of O6-ethylguanine as molecular dosimeter for extrapolation from cells in culture to the mouse. Environ. Health Perspect. 62:163-169. Vousden, K. H., J. L. Bos, C. J. Marshall, and D. H. Phillips. 1986. Mutations activating human c-Ha-rasl protooncogene (HRAS1) induced by chemical carcinogens and depu- rination. Proc. Natl. Acad. Sci. USA 83:1222-1226. Wilson, A. G. E., H.-C. Kung, M. Boroujerdi, and M. W. Anderson. 1981. Inhibition in vivo of the formation of adducts between metabolites of benzo(a)pyrene and DNA by aryl hydrocarbon hydroxylase inducers. Cancer Res. 41:3453-3460. Yang, L. L., V. M. Maher, and J. J. McCormick. 1980. Error-free excision of the cy- totoxic, mutagenic N2-deoxyguanosine DNA adduct formed in human fibroblasts by (+/ - )-7 beta, 8 alpha-dThydroxy-9 alpha, 10 alpha-epoxy-7,8,9,10-tetrahydro- benzo(a)pyrene. Proc. Natl. Acad. Sci. 77:5933-5937.

Resources Available for Simulation in Toxicology: Specialized Computers, Generalized Software, and Communication Networks Daniel B. Menzel, R. L`. Wolpert, ]. R. Boger Ill, and ]. M. Kootsey INTRODUCTION Pharmacokinetics (PK) is a portion of the larger and more general effort of mathematical simulation of physical and biological phenomena. PK simulates certain aspects of the behavior of biological systems, and makes predictions about the future behavior of those systems, by solving systems of algebraic or differential equations. We describe here some of the general aspects of simulation, the com- puter facilities needed to carry out simulations, the specialized resources available nationally for simulations with emphasis on PK, computer lan- guages and approaches used in PK, and some projections for the future development of PK and other biological simulations. Our perspective is that of the toxicologist as a user of these resources. We hope to introduce toxicologists to the field of PK modeling, and to inform those already experienced in PK modeling about some national resources available for their use. G EN ERAL APPROACH ES TO Sl M U LATION Simulation is the representation of physical or biological systems by a set of mathematical relations which approximate reality sufficiently well that predictions from the simulation or model closely approximate obser- vations. The emphasis of this conference was on PK. PK models aim to provide an accurate description of the dose of a drug or chemical reaching 229

230 DANIEL B. MENZEL ET AL. a specific organ and to predict the evolution of that dose with time. This type of PK model (a dosimetry model) is but the beginning of a larger set of models that predict both the dose and the biological outcome, including estimates of risk. Risk estimates fulfill the ultimate objective of toxicology: predicting the health risks to humans (or other potential target organisms) from exposure to toxic chemicals. Most PK models describe the dose of a chemical reaching an entire organ with time. Some subdivide the organ into component parts and describe the chemical reactions within the organ or within individual cells as well as the overall dose of chemical in the organ, whereas others aggregate groups of organs into single compartments for the sake of analy- sis. Regardless of the details, these simulations and PK models all rely on numerical solutions either to partial differential equations or to systems of ordinary differential equations. The problem of numerically solving partial differential equations is usually reduced to that of solving coupled systems of ordinary differential equations by the use of finite-difference methods or, more recently, finite-element methods. Recent progress in numerical analysis has yielded methods that can permit microcomputers to produce useful approximations, even to such complicated systems as the three-dimensional convection of reactive gases in the human lung. COMPUTER LANGUAGES USED IN PHARMACOKINETIC MODELS Any computer language can be used to simulate the distribution ot chemicals within the body or cells with time. The language chosen depends on the computation equipment available, the experimenter's experience and preferences, and the experimenter's objective for the simulation. Early simulation programs were written in general purpose computer languages, but in recent years special languages have been developed for simulation programs to simplify the process of program development and to reduce the amount of specialized knowledge required. Stand-alone models developed on large computers have traditionally been written in the general purpose language FORTRAN. Efficient op- timizing compilers for FORTRAN are available for nearly all computers. There are widely distributed libraries of FORTRAN subroutines for per- forming such common operations as solving algebraic and differential equations, solving linear programming problems, performing elementary optimization, and plotting graphs on video screens or flatbed plotters. Unfortunately, differences among the compilers make it difficult to transport (or "port," in the vernacular) FORTRAN programs and models developed on one machine and operating system to another machine or operating system. In addition, the complete absence in early versions of

Simulation Resources 231 FORTRAN of control structures such as "IF . . . TtIEN . . . ELSE," "DO . . . UNTIL," and "WHILE" loops led to a tradition of writing code which is cryptic and clever and difficult to maintain or upgrade, particularly after the original programmer has left the project. FORTRAN programs require considerable expertise to write and maintain and a major investment in time to develop and debug. More recent languages such as C and Pascal were designed to enable programmers to write even large systems that are easy to maintain. These languages, both derived from the language Algol, have sophisticated con- trol structures that make a program's logical structure and flow of control stand out. This makes it easier for older programs to be understood and maintained by their authors or by other programmers. These modern lan- guages also encourage modular programming, in which each large task (such as simulating the evolution of the concentrations of a chemical toxicant in each of several human organs) is broken down into smaller tasks (such as reading a data base of organ descriptions, building a data structure for the simulation, approximating the evolution of a single organ for 1 s, stepping through the list of organs, stepping through time, and plotting the results), each of which can be broken down into tasks that are smaller still. Each task is performed by a separate subroutine that can be modified and improved as the simulation model evolves, without re- quiring that other subroutines be changed. Experience has shown that large programs are much easier to develop and debug when they employ the principles of modular programming. The languages C and Pascal differ in expected ways, given their origins. Pascal was designed (by the Swiss computer science professor Niklaus Wirth) as a teaching language. It is easy to write well-structured modular programs in Pascal, and very hard not to the language strictly enforces all sorts of rules and gives informative error messages at the first sign of an error. In contrast, the language C was written at Bell Laboratories (by UNIX pioneer Dennis Ritchie) for expert systems programmers. It gives the programmer more power and flexibility than any other high-level language (and nearly as much as assembler language provides). The price paid for this flexibility is that C is harder to learn than Pascal for novice or even for moderately experienced programmers, and programming errors tend to be caught later (when the program is run) rather than earlier (when the program is complied and when errors are easy to correct). Thus, C is preferred for the experienced computer programmer, and Pascal for the novice to intermediate programmer. BASIC, the traditional language of microcomputers, is hardly ever used on larger computers because of its relative lack of control structures and because of its traditional implementation as an interpreted (rather than compiled) language. Interpreters for BASIC are widely available for

232 DAN ~ EL B. M ENZEL ET AL. microcomputers, often with good support for on-screen graphics, and the language is not difficult to master. This makes BASIC appealing for novice programmers and adequate for the smallest simulation projects. C and Pascal compilers are also available for microcomputers, minicomputers, and mainframe computers. The dramatic commercial success of Borland's inexpensive "Turbo Pascal~" for International Business Machines (IBM) personal computers (PCs) and compatible computers has made Pascal the most popular structured language used on microcomputers in general and, therefore, in PK models, whereas the success of UNIX in the minicomputer and engineering workstation world has made C the dominant language on mid-sized computers. Language is hardly ever the limiting factor in PK modeling. SIM ULATION LANGUAGES Simulation languages were developed to provide the constructs and operations required in system simulation. Using such a language, the model developer has only to describe the system to be simulated and does not have to be concerned with writing detailed code for such operations as changing parameters, controlling simulation runs, solving the equations, and printing or plotting the output. Table 1 lists some of the more popular simulation languages; a more detailed listing can be found in the October 1986 issue of the journal Simulation. Because engineers have been using simulation longer than biomedical researchers, most of the available sim- ulation languages are designed for engineering applications and are not ideal for biomedical problems. Many simulation languages, such as CSSL and ACSL, work by trans- lating the simplified model description into a general purpose computer language such as FORTRAN. ACSL is being incorporated into a more TABLE 1 Popular Simulation Languages Language ACSL Description SCoP SIMNON ADSIM A language designed for modeling and evaluating the performance of continuous systems described by time-dependent, nonlinear differential equations. CSSL IV A model development language and translator for solving systems of differential equations and a run-time control language and interpreter. An interactive control program for simulation calculations. Used with a C compiler, SCoP greatly simplifies the construction of a simulation program. Available for both IBM PC and VAX systems. A command-driven interactive program written in FORTRAN for simulation of systems governed by ordinary differential equations and difference equations. A high-level simulation language providing support for the Applied Dynamics International simulation processor, including run-time commands.

Simulation Resources 233 comprehensive program, SIMUSOLV@, which allows even greater flex- ibility and more complex modeling. Elsewhere in this volume, Gary Blau and Richard Reitz discuss PK models that were programmed using SI- MUSOLV. SCoP (Simulation Control Program) is a simulation package developed specifically for biomedical researchers by the National Biomedical Sim- ulation Resource (NBSR) at Duke University. SCoP uses techniques de- veloped for microcomputers such as menus and screen editing to make it easy for both novice and experienced users to interact with a simulation program. Because biological simulation has not yet developed to the stage where models can be constructed from a set of standard modules, SCoP lets the modeler write equations of any type, e.g. linear or nonlinear algebraic, ordinary, or partial differential equations. To build a simulation program, the modeler makes a copy of one of the template files provided and fills in the model equations in the places specified in the file. For clarity, these equations are written in terms of familiar variable names. The SCoP package includes a library of solvers for several kinds of equa- tions, algebraic and time-dependent, as well as numerous other functions to simplify the description of models and experimental setups. The modeler must also build a data base of parameter and variable information using a spreadsheet-like program provided with SCoP. This data base can later be modified interactively through the simulation program. SCoP provides interactive facilities for changing model parameters; selecting, plotting, printing, and saving output; comparing model output with experimental data visually; and adjusting parameters to fit a model to experimental data. SCoP runs on microcomputers such as the IBM PC/XT/AT and on min- icomputers such as the VAX (Digital Equipment Corporation) under the VMS and UNIX operating systems. Thus, a model developed on a microcomputer can be moved up to a larger machine when it gets too large or too slow. EQUIPMENT NEEDS FOR PHARMACOKINETIC MODELING Computational equipment needs vary with the nature of the problem. Complex simulations involving large data sets and multiple equations are best approached using mainframe or minicomputers. Specialized com- puters provide the necessary resources for solution of these complex pro- grams. Conventional multicompartmental PK models can generally be solved on microcomputers, or microcomputers can be used to access larger shared resources. A microcomputer operating as a smart terminal, as well as a remote data processing unit, is the most desirable approach, providing the flexibility of access to shared resources and local computation and data storage for those models not requiring shared resources.

234 DANIEL B. MENZEL ET AL. Specialized computers are available through the Division of Research Resources (National Institutes of Health ENIH]) program in computer resources. Microcomputers equipped with a modest memory of 256 kilobytes are generally adequate for most PK models. The decrease in cost of hard-disk storage makes the availability of 20-30 megabytes of storage accessible to most laboratories. These microcomputers provide adequate computa- tional ability and storage for the majority of PK models, but speed of computation is sacrificed. The development of math coprocessors (such as the 8087/80287 and 68881 chips) and accelerator boards has made PK model solution on microcomputers an achievable objective. Graphic output is generally not a major problem with PK models. Statistical and plotting routines are either incorporated into simulation languages or can be used as general utilities once the PK model has been solved and the results of the prediction obtained. Commercially available packages such as SAS/GRAPH~ Microsoft Chart or LOTUS 1-2-3~ can be used with computer files generated by the PK model program to provide a visual display of the mathematical relationships. A clear graphical rep- resentation of pharmacokinetic relationships is important for conveying model predictions to most toxicologists without extensive training in math- ematics. Public and lay decision makers also understand graphic displays of information more readily than tabular results or equations. Because of the popular need for graphic displays of data, many commercially available programs on microcomputers and minicomputers afford toxicologists the luxury of dispensing with the plot subroutines once required only a few years ago for computer output. Microcomputers or terminals equipped with a modem can easily com- municate with minicomputers and mainframe computers so that specialized computer facilities are within reach of the toxicologist's laboratory. Com- mercially available terminal emulator programs allow experimenters to transport input and output data files from microcomputers up to larger computers and back, so that remote storage and manipulation of PK model output is routine. Programs such as Crosstalk, Smartcom, or PC-Talk, for example, allow the experimenter to assemble on one's own micro- computer large data sets for computing on the larger machines and to store the larger computer's output in one's own microcomputer in return. Graph- ical terminal emulator programs such as PC-Plot III also allow the local microcomputer to display mainframe-generated graphs and charts. CONVENTIONAL APPROACHES TO MODELING The traditional approach to modeling has been the construction of a model which stands alone. Generally, these models have been built to

Simulation Resources 235 describe a complex event and are supported on large general purpose computers (mainframe computers). The regional lung dosimetry model for reactive gases developed by the U.S. Environmental Protection Agency Health Effects Laboratory (Miller et al., 1978, 1985) is an example of such a model. These models are generally written in FORTRAN, C, or Pascal, but rarely in BASIC. Such models can use optimal file formats, data structures, and numerical methods for the needs of a particular model, making them ideal for accommodating large input data sets, complex relationships among the model's state variables, and a need for rapid computation. Mainframe computers running such stand-alone models can provide more memory, speed, and precision than can minicomputers or microcomputers; but most PK models in common use involve few enough equations and are well enough conditioned that the capacity, speed, and precision available on microcomputers are adequate. Especially when mi- crocomputers are used, the user should include diagnostic tests to ensure that rounding errors and other sources of imprecision have not compro- mised the model's prediction. As discussed elsewhere in this volume by Frederick J. Miller and colleagues, sensitivity analysis of the model is essential to ensure that such computational errors do not bias the model results. The reactive gas lung dosimetry model is now available also in a microcomputer version with little loss in precision compared with that supported on mainframe or minicomputers such as the VAX 11 series (F. J. Miller, U.S. Environmental Protection Agency, personal commu- nication). AN EXAMPLE OF A MODULAR PROGRAM The conventional multicompartmental model based on the physiological properties of organs and the chemical and physical properties of chemicals often follows the form shown in Figure 1. This model was developed in our laboratories by Professor R. L. Wolpert and J. R. Boger III. The model program was written in the C language using modular programming techniques. It is a generic model in the sense that it can be used for simulating the time course of chemical concentration in selected organs or parts of organs for different chemicals and for different animal species without reprogramming or recompiling the model program. The model maintains a physiometric data base for each of several animal strains and species, including blood flow rates, organ masses, capillary or plasma volumes, etc., for each animal. For each new chemical the investigator must choose which organs or parts of organs to represent with separate compartments and which to aggregate together. The rates at which the chemical is introduced, elim- inated, or metabolized in each compartment and the extraction ratio re-

236 DANIEL B. MENZEL ET AL. iINHALE: _< E L N o o o S D Or ~ LUNG , LYMPH NODE S _ A SPLEEN ~ R B E L R O L | LIV E R | . | GUT - 1 ,<3 ~ E it' ~ r E S T E it KIDNEY _ ~ t u R N r MU SC L E I | OTHER | FIGURE 1 A generic multicompartmental model for distribution and metabolism of xenobiotic compounds in mammals.

Simulation Resources 237 rating the equilibrium concentrations in the blood or plasma and the tissue of each compartment must also be specified. The model will then plot a color overlay graph showing the time course of chemical concentration in each of a selected number of compartments. Using this model and measured nickel extraction ratios between each organ and plasma, the kinetics of distribution of nickelous ion in the body have been simulated (Francovitch et al., 19861. Extraction ratios were determined in the whole rat by continuous infusion of radiolabeled nickel chloride until a steady-state concentration was attained in blood. A bolus loading dose was used to accelerate an approach to steady-state concen- tra-tions. This model provides a capillary blood volume within each organ within which the tissue and plasma concentrations are at equilibrium following the transit of blood through the organ. Examples of the simulation of a bolus intravenous injection (shown in Figures 2-5) compare favorably with experimental results (also shown in the figures). Lungs, kidneys, and testes are the three organs of the body known to have a particular affinity for nickel. All of the other organs are lumped together to simplify the interpretation of the results, but the data are not lost for individual organs not reported. Such a model is generic in the sense that one simulation model can be used for multiple chemicals and can be simplified according to the bio- z <a ~ 2.5 - - 1 - 2 - - - ~ . o 2 4 TIME (HOURS) 6 FIGURE 2 Simulated amount of nickel in the kidneys following an intravenous injection of 13 ,ug. Experimental results are shown as (O).

238 DAN ~ EL B. M ENZEL ET AL. 200 190 180 1 70 '60 450 1 40 130 1 20 ~ ~ o c, 100 Z go So 70 60 z 50 40 30 20 0 o o i, - - - - - 6 2 4 TIME (HOURS) FIGURE 3 Simulated amount of nickel in the lung following an intravenous injection of 13 leg. Experimental results are shown as ( O ). 900 - BOO - 700 - _% 600- z boo- 1 , 400 - ~ \ 300- 200 - ~ 100 - ~ = o- o 2 4 TIME (HOURS) 6 FIGURE 4 Simulated amount of nickel in the liver following an intravenous injection of 13 lag. Experimental results are shown as ( O ).

Simulation Resources 239 1 ,9 - 1 .8 — 1 .7 - 1 .6 - 1 5- 1.4 - 1.3 - 1.2- 1.1 - 1 - 0.9 - 0.8 0.7 0 6 C] o LO A z 0.5 0.4 0.3 0.2 0 1 - - - I , . , 2 4 Tl ME (HOURS) FIGURE 5 Simulated amount of nickel in the blood following an intravenous injection of 13 log. Experimental results are shown as (O). logical activity of the chemical. An individual organ (such as the kidney or lung) can be represented by several compartments if this is necessary to reflect accurately the transport and metabolism of the chemical being studied. There is little point in reporting and studying concentrations in all parts of all organs for all chemicals when biologically each chemical may affect only one or two target organs. The model can be modified for other species and for different physi- ological data within a species by simply modifying the input data set with a text processor. Allometric relations have been used to construct tables of organ properties within species and between species where direct in- formation is not available in the literature, so that variations in animal size can be accommodated. The development of a standard data set for a 300-g rat has been completed and is shown in Table 2. Models for children, adult human females, pregnant human females, and senescent humans are projected. Similar models for mice and rats of both sexes and different stages of maturation and senescence are under way. AN EXAMPLE SCoP PROGRAM As an example of the application of SCoP to toxicology, a simulation of the effects of sulfite or sulfur dioxide on the covalent reaction of benzoLa~pyrene with DNA or its detoxification by conjugation to unreac-

240 DANIEL B. MENZEL ET AL. TABLE 2 An Example of Physiological Data Needed for a Generic Model: Male 300-g Sprague-Dawley Rat A. Anatomical Compartment Dimensions Capillary Extraction Elimination Compartment Mass (g) Volume (ml) Ratio (ml/g) Rate (liter/in) Lung 1.27 0.126 0.86 0.0 Kidney 2.10 0.013 17.4 0.85 Liver 11.71 0.090 0.31 0.025 Gut 4.60 0.016 0.36 0.0 Spleen 0.54 0.0046 0.36 0.0 Testes 2.73 0.0022 0.31 0.0 Muscle 150.0 0.194 0.18 0.0 Carcass 111.46 0.516 0.33 0.0 Venous blood o.oa 12.90 1.0 0.0 Arterial blood o.oa 4.30 1.0 0.0 B. Blood Flow Rates Among Anatomical Compartments Source Destination Flow Rate (ml/h) 1,241 1,522 947 237 61 1,795 1,241 338 947 237 61 1,795 846 846 175 175 5,640 5,640 Kidney Liver Gut Spleen Testes Carcass Arterial blood Arterial blood Arterial blood Arterial blood Arterial blood Arterial blood Arterial blood Muscle Lung Arterial blood Venous blood Lung Venous blood Venous blood Liver Liver Venous blood Venous blood Kidney Liver Gut Spleen Testes Carcass Muscle Venous blood Venous blood Lung Lung Arterial blood aMost compartments contain both a tissue (mass) and blood (capillary volume) compartment. However, the blood compartments are unusual in that they contain no tissue. live metabolites has been completed (Keller et al., in press). In Figure 6 are shown the chemical reactions involving the production of the putative ultimate carcinogen of benzofa~pyrene, benzo~a~pyrene diol epoxide (BPDE), and its reaction with nuclear DNA to form covalent adducts. 13PDE is conjugated with glutathione (GSH) by glutathione S-transferase

Simulation Resources 241 ~ I So 2- SO2- 3CNU .h AH sp ~ Ow G' SO ~ X a ~ D \m —7 ~: ASH MICROSOMES ~1 C ~ BPOE GET t~.,\~\'\ I\ . GS-BPD ~ 6860f l 1 FIGURE 6 The overall scheme of benzo(a)pyrene chemical reactions involving the metabolism of benzo(a)pyrene to benzo(a)pyrene dial epoxide, and the reaction of benzo(a)pyrene dial epoxide with glutathione and nuclear DNA to form covalent adducts. The detoxification of benzo(a)pyrene diol epoxide by glutathione is inhibited by glutathione S-sulfonate, which is formed by a reaction between sulfite and glutathione disulfide and by added BCNU. to form the unreactive dial conjugate. This pathway effectively removes BPDE from reaction with DNA. The glutathione S-transferase pathway is the rate-limiting step in the detoxification of most polyaromatic hydro- carbon carcinogens. In lung tissue, alternative pathways are limited, mak- ing the lung highly dependent on the functioning of the glutathione pathway. Even so, the lung has limited stores of glutathione. Exposure of cells or the lung to sulfur dioxide or bisulfite results in the formation of glutathione S-sulfonate (GSSO3H) from the naturally occur- ring glutathione disulfide (GSSG). GSSO3H is a competitive inhibitor of glutathione S-transferase (Leung et al., 1985, in press). Exposure of cells to sulfite results in the competitive inhibition of this major detoxification pathway of benzofa~pyrene LB(a)P]. GSSO3H and GSSG are reduced by the same enzyme glutathione re- ductase (Leung et al., 1985~. These two substrates compete with each

242 DAN ~ EL B. M ENZEL ET AL. Other for reduction by glutathione reductase. The Km and Vma,` values for the two substrates differ significantly. Production of GSSG also depends on the availability of reducing equiv- alents within the cell. Reduced nicotinamide adenine dinucleotide phos- phate (NADPH) is the principal source of reducing equivalent, but NADP + may be increased by the flow of reducing equivalent through other path- ways. NADPH is also consumed by other pathways within the cell. At present, it is only possible to estimate the oxidation of GSH to GSSG in general. Experimental values for GSH and GSSG serve as the basis for choosing the initial values in the cell. The reactions described are interdependent, and some compete with each other. Thus, the overall fate of B(a)P through metabolism, detoxi- fication, and reaction with DNA is not easily predictable from a qualitative or intuitive analysis. The rates of these reactions can be described math- ematically, including the efflux of GSSG and GSH from cells by the equations shown in Table 3. The form of the equations suitable for sim- ulation using SCoP are shown in Table 4. The ease of conversion from conventional mathematical form to that needed by the program is a major strength. The names of variables can be made more logical through the use of words descriptive of the variable. The integrated form of the Michaelis-Menten equation is used to de- scribe the enzymatic reactions. To increase the precision of measurement TABLE 3 Equations Describing Benzotaypyrene, Glutathione, and Sulfite Metabolism and DNA Adduct Formation = k3[gssg] [SO3] + 2[Vgssg] + VgssO3h + k5[S03][ h] — k, [ash] — Vex — Vegsh dt -= k,[gsh]l2—k3[gssg][S03] — Vgssg dt FIGS ' = k3[gssg][S03] — V~ssc3h—k20[gssO3h] bEPOXIDE dt bADDUCT dt = = k7[diol] — Vex—k2[epoxide][DNA] - = k2[epoxide][DNA] - k4[adduct] V = - Vgsso3h Vex = Vm gssg [gssg] gss8 Kmgssg + [gssg] VmgssO3h [gSSO3h] KmgssO3h + [gSSO3h] Vmex [ex] Kmex( 1 + [gssO3h]/K;) + [ex] Vegsh = klO(e( k8 t) + k, l (eon t)) VegsSg = kl4(e( kl2 ,) + kI5(e( kl3 i)) · VegSsg + peroxidase

Simulation Resources 243 TABLE 4 Modifications of Mathematical Equations Describing SO2-Benzo (a~pyrene Interaction in a Form Usable by the Simulation Language SCoPa D ash = k3 * gssg * so3 —kl * ash + 2 * V gssg + V_ gsso3h V_ex + kS * so3 * (1/gsh) — V_egsh; D gssg = kl * ash / 2—k3 * gssg *so3 — V gssg— V_egssg + peroxidase; D gsso3h = k3 * gssg * so3 — V gsso3h—k20 * gsso3h; D ex = k7 * diol—V ex—k2 * ex * dna; D adduct = k2 * ex * dna—k4 * adduct; V gssg = V mgssg * gssg I (Kmgssg + gssg); V gsso3h = Vmgsso3h * gsso3h / (Kmgsso3h + gsso3h); V_ex = Vmex * ex / (Kmex * (1 + gsso3h / Ki) + ex); V_egsh = klO * exp( - (k8 * Time)) + kl 1 * exp( - (k9 * Time)); V_egssg = kl4 * exp( - (kl2 * Time )) + klS * exp( - (kl3 * Time)); aSee Table 3 for conventional mathematical representation. Of the components of the GSH pathway, the amounts of GSSG and GSSO3H were increased during the experiment by prior treatment of the cells with the glutathione reductase inhibitor BCNU. By measuring the Km, Ki, and Vmax for all of the enzymatic reactions, the simulation can be solved for the condition with or without BCNU inhibition. In other words, the ab- normal physiological state of inhibition by BCNU can be used to predict the normal physiological results occurring in the absence of BCNU. Actual experimental values were measured in the absence of BCNU to confirm the validity of the assumptions made for BCNU inhibition. Differences in chemical reactivity, such as between syn and anti isomers of BPDE with DNA and glutathione S-transferase substrates, can be ac- counted for also. Experimental values for these rates can be used in the compiled form of the PK model with SCoP and can be modified inter- actively without having to recompile the model. Interaction between the modeler and the model is made much more rapid by this technique. An example of the time course formation of the syn and anti BPDE adducts with DNA over the experiment is shown in Figure 7. The effects of sulfite treatment are shown by the increased B(a)P binding to DNA throughout the course of the experiment. The results agree well with the experimental values found by Leung et al. (in press). Increased B(a)P- DNA adduct formation in the presence of sulfite or sulfur dioxide is mostly due to the inhibition of the glutathione S-transferase pathway by the sulfite metabolite GSSO3H. Using this model and a simple calculation of the intracellular sulfite concentrations that are likely from exposure to sulfur dioxide concentra- tions used in three studies of the cocarcinogenicity of sulfur dioxide with

244 DAN ~ EL B. M ENZEL ET AL. z Ad m Cal - - m Be in AS Cal Be TO O' 10 - O- 0.5 me SUUITE 5.0 me SUUITE (a) - —5 - 30 - A z m 20- - - m at tt' 10- LO a: Z 0- _o o— an n mM SU[FITE - I 1 1 1 1 1 1 1 T i I I I I I · I I 1 5 10 15 20 25 30 35 40 45 50 55 60 MINUTES 0.5 me SUUITE s n mM ';ULFITF 20.0 me SULFITE (b) - I I I I I I I I 1 1 ~ 1 0 5 10 15 20 25 30 35 40 45 50 55 60 MINUTES FIGURE 7 Formation of the anti (a) and syn (b) BPDE adducts with DNA. polyaromatic hydrocarbons (Laskin et al., 1976; Pauluhn et al Pott and Stober, 1983), we can predict that a marked increase in DNA adducts occurred in the study done by Pauluhn et al. (1985), in which there was a significant increase in the number of lung tumors and shortened time to tumor appearance. Only very small increases of 3-4% in DNA adducts were predicted to have occurred in the other two studies. The changes in numbers and time to tumor for these studies was marginal and of questionable statistical significance (Figure 81. ., 1985;

Simulation Resources 245 _ _ 30- o At L' 20- o At a m ~ 10- cn o // / l / I I / I / I i/ / - Laskin ot al. (10 ppm) Pott and Stober (11.5 ppm) Pauluhn et al. (1972ppm) ~ 50 100 150 200 250 300 350 MINUTES FIGURE 8 Predicted increase in amounts of syn-BPDE adduct over control from studies by Laskin et al. (1976), Pott and Stober (1983), and Pauluhn et al. (1985).

246 DAN ~ EL B. M ENZEL ET AL. Our model is still only a crude approximation of the actual flux of BAP through the lung cell and the effects of sulfite or sulfur dioxide on these reactions. The relationship between increased DNA binding and tumor rates is unknown at present. The rate of removal of adducts under these conditions has not been measured. The model overestimates the adduct formation as a consequence. SHARING RESOURCES As stated at the outset of this paper, the simulation of toxicological events, either pharmacokinetics or pharmacodynamics, is but a special case of the more general methodology of simulation of physical, chemical, and biological events. The same mathematical approaches are applicable to all of these fields, as are the same computational resources. Fortunately for toxicologists much of the developmental work has been accomplished in other fields, particularly engineering, and is now available for appli- . . . cation In toxins ogy. The sharing of resources, especially computational resources, provides greater power to pharrnacokinetics and pharmacodynamics than could be justified by the field alone. The NBSR provides a major computation facility specifically designed for biological simulation (see Table 51. Training in simulation is available from the NBSR and other NIH simulation resources at two levels. Introductory classes are designed for researchers with little or no background in computer usage or program- ming. The NBSR presents introductory classes lasting four and a half days several times a year at locations where there is a group of 8 to 15 students and a suitable cluster of microcomputers for teaching. The NBSR intro- ductory class is built around the SCoP software package, and students can take a copy of the software home with them to continue their work after the class is over. Researchers with no prior programming experience can be developing their own simulation programs by the end of the class. Advanced simulation courses are also available or in preparation by TABLE 5 Computation Facilities at the National Biomedical Simulation Resource VAX 11/750 Array Processors (MAP-6430 and Mini-MAP) Special Parallel Coprocessors for Simulation (AD-10 and AD-100)

Simulation Resources 247 NBSR. These courses are for the researcher with some previous experience and cover such topics as the use of special purpose processors, a wider range of simulation languages, and advanced numerical techniques for equation solving. TH E TOXI N CONCEPT To facilitate the exchange of information in the toxicology community, the Toxicology Information Network (TOXIN) has been established within the NBSR. The goals of TOXIN are to improve communication and understanding among toxicologists about mathematical modeling, to pro- vide access to PK and risk assessment models, to develop consensus on PK and risk assessment models through collaborative use of the same model on the same computer, and to provide access to special data bases needed for modeling. A schematic diagram of the relationship between users and TOXIN is shown in Figure 9. Users in the United States can interconnect with TOXIN through the commercial data communications network TYMNET. Users enter TYM- NET through a local telephone number using regular telephone lines. Microcomputers are encouraged as smart terminals for use with TOXIN, but any terminal can be used. Certain terminals have limited graphics, which may prevent the use of some programs. Depending upon the needs of the user, either a TOXIN account or a regular NBSR collaborators account provides access to the facilities of the NBSR, including the spe- cialized computing facilities. Electronic mail for exchange of information, models, and data files is available. PK models contributed to TOXIN are placed in common computer files available to every TOXIN user. The models can then be used to run simulations with the same computer code for use on the same computer by different investigators. Problems, improvements, enhancements, and comments on the models can be made by users and sent to the contributor of the model or to the general community of TOXIN users through the electronic mail system. In this manner a consensus can be built about the critical aspects of PK models. The selection of physiological parameters for PK models is a particularly critical area in which comparisons of different values may prove especially useful. As better physiological val- ues are amassed, they can be stored in data tables in a form useful for modeling. Alternative models can be compared and potentially integrated to form second- and third-generation models. As more PK models are developed and made available to TOXIN users, an archive of PK and pharmacodynamic models can be accumulated. Historical reference will then be possible. At present such an archive does

248 DANIEL B. MENZEL ET AL. National & International Users , 1 | TYMNET 1 . ~ NBSR HOST I TOXIN .1 SPECIALIZED SYSTEMS FIGURE 9 Schematic diagram of the relationship between users and TOXIN. not exist, and original computer codes are inaccessible to investigators beyond the originator. If PK models are to have a major impact on reg- ulatory decision making, the establishment of an archive and the requisite data for risk estimation is essential for the future review and improvement of regulations. Selection among different PK models is now difficult to make because direct comparisons are generally not possible. Storage of multiple PK models in a form that can be easily used by investigators will provide a means of systematic testing and further validation of PK modeling ap- proaches. Estimation of errors within different PK models can be made more easily, and the methodology of estimating error in PK models (as discussed by Gary E. Blau and W. Brock Neely, this volume) is made more feasible. TOXIN also provides access to commercial programs that are expensive and difficult to justify for a single laboratory. At present, ACSL, CSSL, SIMNON, ADSIM, SCoP, and SCoPFIT are available as simulation lan- guages on NBSR and are accessible through TOXIN. Should user demand

Simulation Resources 249 warrant other commerical software, the languages supported by NBSR can be expanded. FUTURE TRENDS The rapid development of simulation languages should make simula- tions of all kinds much more accessible to nonprogrammers. These lan- guages have gone through several revisions and are improving rapidly. Reductions in costs are likely as the programs become more widely used. Versions for microcomputers will provide even greater access. Even those programs not supported by microcomputers should become more available through the use of shared resources and the telecommunications network. Shared resources provide the opportunity to accumulate libraries of PK models that can be used by multiple users. Consensus can be built and improvements in shared PK models can be accelerated through comparison of predictions for a variety of compounds. Development of special data bases for physiological parameters can be achieved. Much of our knowledge of organ blood volumes, flow rates, and gross anatomy are gained from tracer techniques. Pharmacokinetic data provide similar data, and the concordance of several experimenters' data increases our confidence in these basic physiological values. Errors are likely to be detected when reference values fail to predict experimental results. Data gaps will be identified and filled as physiologically based PK models become more widely used. Computational devices are evolving rapidly. The drastically reduced costs of data storage and the ease of transmission of masses of data by electronic means will result in more complex models useful for experi- mental design and, with hope, long-range planning. Fidelity with human experience should increase as PK models are integrated with pharmaco- dynamic models, making risk assessment more reliable. Distributed data storage will also improve data utilization and decrease repetition of ex- periments, lowering costs and animal usage and increasing generalization among classes of chemicals. Digital analog devices and their high-level simulation languages should speed computation to real time rates. Digital equipment reduces the uncertainties and difficulties of conventional wired analog or hybrid computers. CONCLUSIONS In summary, there are no major equipment limitations for most PK models in the United States. The specialized computer facilities such as the NBSR at Duke University and other components of the NIH simulation

250 DAN ~ EL B. M ENZEL ET AL. resources provide more than ample power for computation. Most PK models can be solved on microcomputers using commercial or public domain software. The rapid development of simulation software provides an opportunity for toxicologists to use the software firsthand to simulate their experiments, to assist in design of experiments, and to extrapolate from animal experiments to humans. ACKNOWLEDGM ENTS This work was supported in part by grants from NIH ES01859, ES02916, CA14236, and RP-940-5 from the National Institutes of Health and by a grant from the Electric Power Research Institute. We wish to thank Drs. K. H. Leung and D. A. Keller for sharing their data and models on sulfur dioxide and benzofa~pyrene metabolism and Drs. R. J. Francovitch and C. R. Shoaf for sharing their data arid model on nickel metabolism. Messrs. M. I. Tayyeb and J. Sandy and Ms. K. Wilkinson provided excellent technical assistance. REFERENCES Francovitch, R. J., C. R. Shoaf, M. I. Tayyeb, and D. B. Menzel. 1986. Modeling nickel dosimetry from kinetic measurements in rats. The Toxicologist 6:528. Keller, D. A., K. H. Leung, and D. B. Menzel. In press. Glutathione, sulfite, and benzo(a)pyrene interactions: A mathematical model. The Toxicologist. Laskin, S., M. Kuschner, A. Sellakumar, and G. V. Katz. 1976. Combined carcinogen- irritant animal inhalation studies. Pp. 190-213 in Air Pollution and the Lung, E. F. Aharonson, A. Ben-David, and M. A. Klingberg, eds. New York: John Wiley & Sons. Leung, K. H., G. B. Post, and D. B. Menzel. 1985. Glutathione S-Sulfonate, a sulfur dioxide metabolite, as a competitive inhibitor of glutathione S-transferase, and its re- duction by glutathione reductase. Toxicol. Appl. Pharmacol. 77:388-394. Leung, K. H., D. A. Keller, and D. B. Menzel. In press. Glutathion S-sulfonate inhibition of glutathione conjugation with benzo(a)pyrene epoxides. Toxicol. Appl. Pharmacol. Miller, F. J., D. B. Menzel, and D. L. Coffin. 1978. Similarity between man and laboratory animals in regional pulmonary deposition of ozone. Environ. Res. 17:84-101. Miller, F. J., J. H. Overton, Jr., R. H. Jaskot, and D. B. Menzel. 1985. A model of the regional uptake of gaseous pollutants in the lung. I. The sensitivity of the uptake of ozone in the human lung to lower respiratory tract secretions and exercise. Toxicol. Appl. Pharmacol. 79:11-27. Pauluhn, J., J. Thyssen, J. Althoff, G. Kimmerle, and U. Mohr. 1985. Long-term inhalation study with benzo(a)pyrene and SO2 in Syrian golden hamsters. Exp. Pathol. 28:31. Pott, F., and W. Stober. 1983. Carcinogenicity of airborne combustion products observed in subcutaneous tissue and lungs of laboratory rodents. Environ. Health Perspect. 47:293- 303.