Below are the first 10 and last 10 pages of uncorrected machineread text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.
Intended to provide our own search engines and external engines with highly rich, chapterrepresentative searchable text on the opening pages of each chapter.
Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.
OCR for page 183
PART IV
Uncertainties: Integration
with Risk Assessment and
Resou roes
OCR for page 183
OCR for page 183
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 concentrationtime 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
OCR for page 183
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 modelbuilding 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
OCR for page 183
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 modelbuilding 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
OCR for page 183
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
OCR for page 183
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 twobytwo 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
OCR for page 183
190 GARY E. BLAU AND W. BROCK NEELY
magnitude (Dunker, 19841. In addition, the stateoftheart 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., twodimensional 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
OCR for page 183
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 modelbuilding 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
OCR for page 183
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 LineweaverBurk 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 firstorder 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)
OCR for page 183
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~4labeled 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
trationTime 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
OCR for page 183
240 DANIEL B. MENZEL ET AL.
TABLE 2 An Example of Physiological Data Needed for a Generic
Model: Male 300g SpragueDawley 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 Stransferase
OCR for page 183
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\ .
GSBPD
~ 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 Ssulfonate, 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 Stransferase pathway
is the ratelimiting 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 Ssulfonate (GSSO3H) from the naturally occur
ring glutathione disulfide (GSSG). GSSO3H is a competitive inhibitor of
glutathione Stransferase (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
OCR for page 183
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 MichaelisMenten 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
OCR for page 183
Simulation Resources 243
TABLE 4 Modifications of Mathematical Equations Describing SO2Benzo
(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 Stransferase 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 Stransferase 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
OCR for page 183
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 34% 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;
OCR for page 183
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 synBPDE adduct over control from studies by
Laskin et al. (1976), Pott and Stober (1983), and Pauluhn et al. (1985).
OCR for page 183
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
(MAP6430 and MiniMAP)
Special Parallel Coprocessors for Simulation
(AD10 and AD100)
OCR for page 183
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 thirdgeneration 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
OCR for page 183
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
OCR for page 183
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, longrange 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 highlevel
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
OCR for page 183
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 RP9405 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. 190213 in Air Pollution and the Lung, E. F.
Aharonson, A. BenDavid, and M. A. Klingberg, eds. New York: John Wiley & Sons.
Leung, K. H., G. B. Post, and D. B. Menzel. 1985. Glutathione SSulfonate, a sulfur
dioxide metabolite, as a competitive inhibitor of glutathione Stransferase, and its re
duction by glutathione reductase. Toxicol. Appl. Pharmacol. 77:388394.
Leung, K. H., D. A. Keller, and D. B. Menzel. In press. Glutathion Ssulfonate 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:84101.
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:1127.
Pauluhn, J., J. Thyssen, J. Althoff, G. Kimmerle, and U. Mohr. 1985. Longterm 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.