- Email: [email protected]
A note on the use of a stochastic mammillary compartmental model as an environmental safety model
BULLETIN
OF
MATHEMATICAL VOLUME
BIOLOGY
38, 1970
A NOTE ON THE USE OF A STOCHASTIC MAMMILLARY COMPARTMENTAL MODEL AS AN ENVIRONMENTAL SAFETY MODEL
q J. H. MATIS Institute
of Statistics,
Texas A & M University
College Station,
Texas, 77840, U.S.A. R. L. KODELL Kational Center for Toxicological Research, Jefferson, Ark., 72079, U.S.A.
M. CARDENAS Department of Experiment.al Statistics, New Mexico State University, Las Cruces, N.M., 88003, U.S.A.
Environmental safety testing typically requires procedures for extrapolating from the relatively high experimental to the very low use doses of potentially harmful substances. In the present paper, a stochastic mammillary compartmental model for environmental safety testing is proposed and extrapolation procedures based on its dose-response relationship are developed. The proposed model is a direct generalization of one of the basic safety models, the one-hit model, in that a harmful reaction is assumed to occur if at any time any of the peripheral compartments attains a specified threshold of particles. consideration of a closed model yields an upper bound on the probability of attaining a certain threshold level, thus providing a conservative procedure for extrapolating to a low dose, while a lower bound obtained from a related open model provides a useful monitoring device as t.o the sharpness of the upper bound. The extrapolation procedure is illustrated with simulated data and approximations for initial values are developed.
The determination of low-risk dose levels of potentially 1. Introduction. harmful substances in the food, drug and chemical industries, with particular has been the subject of considerable application to suspected carcinogens, 467
468
J. H. ~MA’I!IS, R. L. KODELL
X3-D 11. CARDESAS
recent interest and .controversy (see e.g., Symposium of Statistics and the Environment in J. Wu&. AC&?. Sci. (1974) and Scielzce (Oct. 1974), 18, 239-245). Most of the ~ontrove~y centers about tec~~~iq~~es of extrapolatioil both from the relatively high experimental to the very low use doses and from “mouse-to-man”. Experiments are conducted with high dosage levels and, in the absence of a known functional relationship between dose and response in the very low dose range, “conservative” extrapolation procedures have been developed as outlined by Mantel and Bryan (1916) and an FDA Advisory Committee Report (1971). The primary purpose of this paper is to propose a stoch~tic mamm~la~ compa~mental model for environmental safety testing and to develop extrapolation procedures based on its dose-response relationship. The paper also explores certain properties of the model but does not consider the details of statistical parameter estimation. The references above outline the three principal safety models used for doseresponse extrapolation for a single pulse dose; the probit, the logit, and the one-hit models. The first two are justified primarily by their adequate approximation of many observed dose-response relationships in the accelerated dose range. However there seems to be little or no evidence, either empirical or theoretical, of the validity of these models in the very low dose range. On the other hand, the one-hit model does represent an attempt to construct a rough causal model for low doses and hence many researchers studying radiation as a carcinogenic process favor the one-hit model of the process while others theorize “that it takes some minimum number of molecules, not a single molecule, to produce an effect” (Schneiderman (1974) p. 72)). The theoretical analysis by Arley (1961) and the analysis of virus data by Bryan (1961) are just two of many studies supporting the one-hit and the multiple-hit models in the case of carcinogenesis. Yet in their present form, extrapolation based on these models leads to extremely conservative results and hence they are not widely used. The ~ompa~mental safety model herein proposed and also previously investigated by Matis, Cardenas and Kodell (1974) generalizes the theoretical appealing one-hit model to a multiple-hit model within the broader framework of compartmental analysis. Compartmental analysis is a technique which has broad acceptance in the biomedical sciences in general and in physiological modelling in particular (see e.g., Atkins (1969), Jscquez (1972), Rescigno and Segre (1966), or Sheppard (1962)). A compartmental model typically does not purport to be a detailed causal model of a phenomenon but instead represents a rough and readily eon~ptualized framework su~ciently rich in parameter to provide an adequate approximation to many physiological systems. Moreover the recent addition of stochastic variation to the compartmental model (Matis and Hartley (1971)) provides a rationale for making valid inferences about a complex phenomenon from the relatively simple compartmental model
USE OF A STOCHASTIC MAMMILLARY
COMPARTMENTAL MODEL
469
inasmuch as the need for the intricate detail of the phenomenon is supplanted by probabilistic variation. The theoretical justification of the one-hit and multi-hit models combined with the general acceptance and versatility of compartmental models make a strong case for the present safety model. Although most applications of the cited safety models deal directly with carcinogenic safety testing, the basic principles apply to a wide variety of other environmental and ecological applications. It should also be pointed out that there exists another class of safety models based on a time to response (tumor) relationship (see e.g., Armitage and Doll (1961), Peto and Lee (1973), and Hartley and Sielken (1975)) which is not considered in this paper. 2. ilfodel I-An Upen M~~nrn~ll~T~ Compartment 3fodel. We assume the mammillary compartmental model illustrated in Figure 1 and designated as Model I. Let Xi(t), for i = 1, , . . , k, denote the number of particles in each of
h. a
1
Figure 1.
K2=
Ia x,(t)i 1”:
r5 a
.
l
0
x,(t)
A simple mammiIl&y comp~rtmen~l compartments
‘c20
i
model w&h k peripheral
the peripheral compartments at time t and Xk+l(t) the number remaining in the central compartment. It is also assumed that the central compartment is pulse lnhelled with a known number, say N, of independent particles at time 0, \vhereupon Xk+i(O) = M and X,(O) = 0, i = 1, , . ., k. The parameters a, ~1, and 7~sdetermine the fixed transition intensity c~~cients according to the following definitions (for small At): Prob. Prob. Prob.
{a, given particle in camp. k + 1 at time t transfers into camp. i in the interval i = 1, . . ., k. (t, t+At)) = a&+0(&), {a given particle in camp. k+ 1 at time t leaves the system in the interval (t, z+At)f = x,aA.1f O(At), {a given particle in camp. i at time t leaves the system in the intervai (t, t+ At)] = K2aAt+ OfAt), i = 1, . . ., k.
These assumptions on t,he transfer p~obab~l~t~es are analogous to the exact transfer proportions specified in t.he more restrictive deterministio model. ,The model is termed “open” in the peripheral compartments for ~2 > 0.
470
J. H. MATIS, R. L. KODELL AND JI. CARDESAS
In a physiological safety model, the central compartment may be envisioned to be t,he blood plasma and the k peripheral compartments are either cells or groups of cells. This conceptualization is suggested by the Jacquez et al. (1960) model of the distribution of drugs in the body. A reaction (e.g. carcinogenesis) may be assumed to occur in the model if at any time any of the k peripheral compartments attains a specified threshold size, say of s particles. The notion of both fixed and stochastic thresholds has been widely used in carcinogenesis, and the above reaction would be characterized as an s-hit, single stage mechanism by Neyman and Scott (1966). Finally, the openness property allows repair within the peripheral compartments. As an example of an environmental application, a cotton field represented by the central compartment and the k adjacent fields by the peripheral compartments. An adjacent field becomes very attractive to and hence also infested with boll weevils if ever a threshold number of s such insects chance to be present at any time. Whatever the application, Model I undoubtedly represents a gross oversimplification of the phenomenon and many generalizations are possible, some of which will be discussed subsequently. Nevertheless Model I as presently constituted has the great advantage of parameter parsimony and yet at the same time provides much richer models for extrapolation than those presently available. It was shown in Matis, Cardenas and Kodell (1974) that the multivariate probability distribution of the random vector Xi(t), X,(t), . . . , Xk+l(t), X,+,(t), where k + 2 denotes the exterior compartment, is a (k + 2)-variate multinomial with parameters N, pi(t), . . . , pi(t), pz(t), pa(t) where pi(t) = (k + ~1 - Ks)-l{exp ( - Kzat) - exp ( - a(k + p2(t) =
exp (- a(k+ K$),
~30)
1 --h-W
=
KI)~)),
and
-pzV).
Letting Xi denote the maximum which camp. i ever attains over time, and Y’ the maximum of Xi, Xk , . . . , XL, it follows that Prob [Y’ 2 s] denotes the probability of a reaction. The analytical expression for Prob [Y’ 2 s] is intractable but it is shown that a lower bound is Prob [Y’ 2 s] 2 1 - {Prob [Xi < s]>~ 2 1-
8-l
c
{ f-0
jy 0 i
where p = (K~)-~c-~‘(~-~) with c = (k + K~)/KZ . Although the above is useful in its own right, a safety model requires an upper bound on Prob [Y’ 2 s] for use in subsequent extrapolation. We propose the procedure below for generating this upper bound.
USE OF A STOCH.4STIC
MAMMILLART
COMPARThfENTrlL
MODEL
471
3. Model II-A Closed Mammillury Compartmental Model. The closed mammiilary comp~tmental model, where ~2 = 0 and consequently where each peripheral compartment is a ‘%rap” from whence particles cannot escape is designated Model II. Lettin g Tg(t) denote the number of particles in compartment i, it follows that Cfz,s T,(t) = N. In the safety model context, a harmful reaction would be assumed to occur if any peripheral compartment reaches a threshold of s cumulative part,icles. Model II, in addition to providing an alternative safety model, also yields an upper bound for Model I as follows. Let T; denote lim T,(t) over time. Since the particles are independent and the eompart~nents mutually exclusive, it follows that the multivariate probability distribution of !l’; , Ti, . . . , l’L+z is also a (k + 2)-variate multinomial with parameters N, pr, ~1, . . . , pi, 0, p3; where pl = (k+ &-I and 13s = ~i(k+ ~$1. Let S’ be the maximum of T; , Th , . . . , Ti, i.e., the maximum peripheral compartment. One may then conclude from t,he Bonferroni inequality (Feller (1968)) that Prob[S’
fe,
h-2,8).
472
J. H. MATIS,
R. L. KODELL
AXD
M. CARDENAS
The entries in Table I illustrate numerically the performance of the bounds, fa and f4, as approximations to the true probabilities, fi and for specified thresholds, 8. Values offi and fz were generated by computer simulations with a = 1, k = 10, rcl = 10, and ~2 = 1 or 0 for Models I and II, respectively. These empirical values offi and f2 are the result of 1000 realizations of the given models for each value of N. The values of f3 and f4 were determined from expressions (1) and (2). For a physiological safety model, the application mentioned above becomes inverted. That is, for a specified, typically very low risk of p*, the low-dose
fz,
TABLE I Numerical Illustration of Probability Bounds N
40
60
80
100
S
f2
fs
f4
fl
5 6 7 8
0.2392 0.0648 0.0139 0.0025
0.341 0.087 0.019 0.003
0.413 0.124 0.030 0.006
0.4803 0.1388 0.0339 0.0071
6 7 8 9
0.3527 0.1297 0.0383 0.0097
0.496 0.204 0.057 0.012
0.612 0.280 0.094 0.027
0.7872 0.2970 0.0979 0.0285
8 10
0.1908 0.0688 0.0214
0.282 0.099 0.034
0.416 0.163 0.059
0.4669 0.1840 0.0654
9 10 11
0.2433 0.1008 0.0363
0.356 0.155 0.058
0.506 0.254 0.111
0.6309 0.2819 0.1147
9
“safe” level of Model I, say nz , or of Model II, say nf , are required, but they are available only through extensive simulation. In their stead, a “conservative” dose nl,could be derived from the upper bound function, f3, and the lower bound function, f4, could be used to monitor the sharpness of the upper bound, since f3(nz)-f4(ng) represents the maximum possible deviation of f3 from the true probability, be it fi or f2. If, on the other hand, the application required a “liberal” dose, the roles of f3 and f4 would be reversed, and a liberal dose, nz, could be derived. Since, clearly, nj s ni 2 n: 2 nz, the difference nz - ng would also provide a convenient monitoring device, whatever the application. The interrelationships of fi, f2, f3 and f4 as functions of iV are represented graphically in Figure 2. Input parameters were those used to generate Table I,
USE OF A STOCHASTIC MAMMILLARY COMPARTMENTAL MODEL
473
the threshold level being s = 8. There appears to be good convergence in the low-dose extrapolation region of the curves, as will be borne out by the following application. Suppose that in an environmental problem, values of s = 8, k = 10, ~1 = 10 are known parameters relating to a certain chemical agent for which a low-dose “safe” level is desired. Assume further that a risk level of 1 in 10,000 haa been ascertained as acceptable. Thus, regardless of which of the two compartmental
DOSE N
Figure 2. The probability of exceeding the threshold under Model I (fi) and under Model II (fi) with upper and lower bounds (A andfa) as a function of dose size N. The parameters for this example were 8 = 8, k = 10, a = 1, KI = 10 and Kz = i or 0
models applies, to determine the value of nl corresponding to a risk level of p* = l/10,000, one solves f&s)
= 10{ 1- ,-,(y)(-L)” i
(LY)“-‘1 = 0.0001
for N = ni . Repeated evaluation offs(N) on a computer for a range of N-values yields f3(23) = 0.00009713 and f3(24) = 0.00013928. A conservative choice of dose for the given risk !evel would then be ng = 23. A check on the precision of
474
J. H. MATIS,
R. L. KODELL
&SD
JI. CARDESAS
the probability bound leading to this choice of close is made by computing f4(23) = 0.00003042. An appreciation of the excellent capability of the extrapolat8ion procedure can be gained if the acceptable risk level is decreased drastically. Consider now an acceptable risk level of only 1 in 100,000,000, a risk level which is common in the safety testing of potential carcinogens. To determine an acceptable dose, ng ,one solves fs(N)=l
$ I-
,_,( ; X&.~(~)“-“)
=0.00000001
for N =nz. As before, Orialand error yieldsfs(9) =3.359x 10Tgandfs(10) = 1.605 x lo-*. Again, a conservative choice for an acceptable dose is nz = 9. The corresponding lower bound for a close of AY= 9 and a threshold level of 8 = 8 is f4(9) = 9.581 x 10-l’, thus indicating the precision of the upper bound on risk. 5. On the Implementation
of the safety Node1
The previous trial and error procedure involved in finding the safe dose, nz, by solving for N in p* =fe(N) can be simplified somewhat by using an initial estimate of 1~: based on an approximation to fa. Let a new constant E.be defined as (a) Initial
approximation
of nz.
#I = l-
p*/k =
s-1 iv 1 . ti=o0 a
pf(1 -p)N-“.
Inasmuch as p is small, one can use the Poisson approximation to the binomial to obtain R N CfZi exp (-Np)(N p )i/‘l z ., w h’lc h CZIAbe rewritten as the cumulative ~2 distribution (see e.g. Beyer (1963)) as A=
s
O”exp
( - x2/w2/w-1
dX2
2r(s)
~NP
Hence 2Np is the upper 1 percentage point of a ~2 distribution with parameter 29, usually denoted x&r+ and therefore nz 2: xzs, 142~. To illustrate the approximation, assume again that s = 8, k - 10,~1 = 10 and p* = 0.0001, from whence p = 0.05, 1 - i, = 0.00001 and n3 =
(x:3,
o.oooolW~.
Using Table II by Lentner (1975) for the extreme lower percentage points of the ~2 distribution, one has ns = 1.9931/0.1 N 20. With this initial estimate, the procedure would rapidly converge on the safe dose nz = 23. In the same
_
_
0.05233104 0.13020332 0.25604775 0.43368150 0.(16383980 0.94544407
0.~000000 0.00002828 0.00168722 0.01401817
0.~~~1
0.08315227 0.19195597 0.35804763 0.58303141 0.86606411 1.20426967
0.00000000 0.00008944 0.00303589 0.02495551
0.00000000~ -
0.13232794 0.28359924 0.50198774 0.78641547 1.13355429 1.53921731
0.00000002 0.00028286 0.00783741 0.04446447
0.000~00~
_F
~I__-
--.-
._
_
2-‘qI-(v/2)1-1
0.22110412 0.42033933 0.70644114 1.06496937 1.49021874 1.97637282
0.00000020 0.00089450 0.01690433 0.07934626
0.0~~01
t Courtesy of Marvin Lentnor, Dept. of Exp. Stat., New Me?tioo &da Univ.
10 12 14 I6 18 20
2 4 G 8
28
II
Percentage Points of ~2 Distribution* F[x~/v] =
TABLE
0.25535777 0.54516954 0.93959477 1.42018445 1.99310094 2.F3012172 3.32867533
0.14195478 0.33812600 O.F2GO(i29G 0.99903807 1.45069332 1.97113791 2.55363768
0.00001 ___
2/2)aY2-1
0.0~002000 0.00895763 0.07907433
-~
(-
0.000002~~0 0.00282970 0.03G50857
0.000001
j”o* exp
ckc
0.40359379 0.88892030 1.42749511 2.06079786 2.77393852 3.55516529 4.39516272
0.00020001 0.02841848 0.17235211
O*OOOl
476
J. H. MATIS,
R. L. KODELL
AND
Bf. CARDEXAS
problem with p* = 0.00000001, it follows that ns = Cxz, o.ooooOOoOl)/O.l = 0.5330/0.1 = 6, as compared to ~2.:= 9. It can be shown that an approximation to 1~: is also -C/C- 1 . 124= x&+JZ~ with 1 = (I- p*) ’ and p = (tc#1 c (b) linknown parameters k, ~1, ~2, 8. In many environmental applications the parameters have a well-defined structural interpretation and they are either known or measurable by suitable ex~rimentation. In physiological applications, however, the parameters may not correspond as closely to recognizable physical characteristics and hence it is less likely that the parameters can be specified initially. As with other safety models, the case of unknown parameters necessitates testing the response at multiple dose levels, say where ~lt exceeds the number of u~no~ parameters. Xl, N2, . . . , N, Equation (I) or (2) may then be fitted to the dose-quanta1 response data and the parameters thereby estimated. Although the estimation procedures raise some auxiliary statistical questions, the same general extrapolation principles hold. The analytical study and simulation of the statistical procedures are a topic of current research. (c) Generalizations. It is not difficult to postulate a multitude of realistic generalizations of this basic mammillary model. One immediate generalization would arise by allowing unequal transfer coefficients both into and out of the ~ompa~ments, say olg and 66, i = 1, 2, . . . , k, respectively, as considered previously by Matis, Cardenas, Kodell(1974). Another natural extension would be to incorporate a stochastic threshold where s would have a probability distribution and one could also plausibly argue that the transition intensity coefficients must almost surely be time dependent and/or that N is a random variable with a known probability distribution (see e.g. Cardenas and Matis f 1974)). Although involved expressions could be written for each of these generalizations, it is likely that the present model is very robust against such alternatives. However as data in the low-risk dose range becomes available, the validity of the basic model will be tested and its structure refined to fit particular applications. 6. Discussion. The incorporation of the multiple-hit reaction mechanism within a compartmental structure suggests a broad class of safety models of which the present safety model using a particular mammillary system is only one special case. Besides their conceptual appeal, those models offer a family of curves for extra~lation from the high-risk experimental doses to the low-risk use doses. It may also be possible to extrapolate from “mouse-to-man” with these models by correlating a set of species specific parameter estimates with
USE OF A STOCHASTIC MAMMILLARY COMPARTMENTAL MODEL
477
some well-defined physiological characteristics, thereby extrapolating from a hierarchy of animal species to man. The referees of our previous paper, R. Elashoff, A. Rescigno and C. W. Sheppard, made several suggestions for its extension, some of which are incorporated in this manuscript. We acknowledge and appreciate their comments on the previous and the present paper. The investigation was supported in part by a Public Health Service Research Career Development Award 1 K04 BM00033-01 from the National Institute of General Medical Sciences to the first author.
LITERATURE Arley, N. 1961. “Theoretical Analysis of Carcinogenesis.” Proo. Fourth Berkdey Symp. Math. Stat. Prob., 4, 1-18. Arm&age, P. and R. Doll. 1961. “Stochastic Models for Careinogenesis.” Proc. Fourth Berk~y Symp. Math. Stat. Preb., 4, 19-38. Atkins, C. L. 1961. ~~~~t~orn~~n~ent Models in Biological Sy8teme. London: Methuen. Beyer, W. H. (ed.) 1968. Handbook of Table8 for Probability and Stati8tie8, 2nd Ed. Cleveland: Chemical Rubber Co. Bryan, W. R. 1961. “Virus Carcinogenesis.” Proc. Fourth Berkeley Symp. M&h. Stat. Prob., 4, 123-152. Cardenas, M. and J. H. Matis. 1974. “On the Stoch&ic Theory of Compartments: Solution for ~-Compa~ment Systems with IrreversibleTime-De~ndent TransitionProbabilities.” Bull. Math. Biol., 36, 489-504. FDA Advisory Committee on Protocols for Safety Evaluation. 1971. “Panel on Carcinogenesis Report on Cancer Testing in the Safety Evaluation of Food Additives and Pesticides,” Toxic. Appl. Pharmac., 20, 419-438. Feller, W. 1968. An Intmduction to Probability Th,ewy and It8 AppEieatiaow, Vol. 1, 3rd Ed. New York: Wiley. Hartley, H. 0. and R. L. Sielken. 1975. “A Parametric Model for ‘Safety’ Testing of Carcinogenic Agents.” FDA Technical Report No. 1, Institute of Statistics, Texas A 6 M University. Jacques, J. A. 1972. Compartmental Andysia in Biology and Medicine. Amsterdam:
Elsevier .
R. Bellman and R. Kalaba. 1960. “The Distribution of a Drug in the Body.” B&E. Xath. Biop~ys., 22, 309-322. Lentner, Y. 1975. Percentage Points of x2 Distribution. Personal Communication. JIant,el, N. and W. R. Bryan. 1961. “Safety Testing of CarcinogenicAgents.” J. Natn. CancerInst., 27, 321-346. Matis, J. H. and H. 0. Hartley. 1971. “Stochastic Compartmental Analysis : Model and Least, Squares Estimation from Time Series Data.” Biometric& 27, 27-102. &f. Cardenasa.ndR. L. Kodell. 1974. “On the Probability of Reaching a Threshold -7 in a Stochastic ~~a~ill~ System. ” Bull. Hath. Bid., 36, 445-454. Neyman, J. and E. Scot,t, 1966. “Statistical Aspects of the Problem of Carcinogenesis.” Proc. Fijth Berkeley Symp. Math. Stat. Prob., 4, 745-776. Feto, R. and F. Lee. 1973. “Weibull Distributions for Continuous Carcinogenesis Experiments.” Biometrics. 29, 457-476. Rescigno, A. and C. Segre. 1966 l?r~cgand Tracer Kinetics. Waltham: Blaisdell.
-,
478
J. H. MATIS, R. L. KODELL AND M. CAFtDENAS
Schneiderman, M. 1974. “Symposium Acud. Science. Sheppard, C. W. 1962. Baeic Primiples
on Statistics
and the Environment.”
J. Wash.
of the Tracer Uethod. New York: Wiley. RECEIVED 7-28-75 REVISED l-8-76
OF
MATHEMATICAL VOLUME
BIOLOGY
38, 1970
A NOTE ON THE USE OF A STOCHASTIC MAMMILLARY COMPARTMENTAL MODEL AS AN ENVIRONMENTAL SAFETY MODEL
q J. H. MATIS Institute
of Statistics,
Texas A & M University
College Station,
Texas, 77840, U.S.A. R. L. KODELL Kational Center for Toxicological Research, Jefferson, Ark., 72079, U.S.A.
M. CARDENAS Department of Experiment.al Statistics, New Mexico State University, Las Cruces, N.M., 88003, U.S.A.
Environmental safety testing typically requires procedures for extrapolating from the relatively high experimental to the very low use doses of potentially harmful substances. In the present paper, a stochastic mammillary compartmental model for environmental safety testing is proposed and extrapolation procedures based on its dose-response relationship are developed. The proposed model is a direct generalization of one of the basic safety models, the one-hit model, in that a harmful reaction is assumed to occur if at any time any of the peripheral compartments attains a specified threshold of particles. consideration of a closed model yields an upper bound on the probability of attaining a certain threshold level, thus providing a conservative procedure for extrapolating to a low dose, while a lower bound obtained from a related open model provides a useful monitoring device as t.o the sharpness of the upper bound. The extrapolation procedure is illustrated with simulated data and approximations for initial values are developed.
The determination of low-risk dose levels of potentially 1. Introduction. harmful substances in the food, drug and chemical industries, with particular has been the subject of considerable application to suspected carcinogens, 467
468
J. H. ~MA’I!IS, R. L. KODELL
X3-D 11. CARDESAS
recent interest and .controversy (see e.g., Symposium of Statistics and the Environment in J. Wu&. AC&?. Sci. (1974) and Scielzce (Oct. 1974), 18, 239-245). Most of the ~ontrove~y centers about tec~~~iq~~es of extrapolatioil both from the relatively high experimental to the very low use doses and from “mouse-to-man”. Experiments are conducted with high dosage levels and, in the absence of a known functional relationship between dose and response in the very low dose range, “conservative” extrapolation procedures have been developed as outlined by Mantel and Bryan (1916) and an FDA Advisory Committee Report (1971). The primary purpose of this paper is to propose a stoch~tic mamm~la~ compa~mental model for environmental safety testing and to develop extrapolation procedures based on its dose-response relationship. The paper also explores certain properties of the model but does not consider the details of statistical parameter estimation. The references above outline the three principal safety models used for doseresponse extrapolation for a single pulse dose; the probit, the logit, and the one-hit models. The first two are justified primarily by their adequate approximation of many observed dose-response relationships in the accelerated dose range. However there seems to be little or no evidence, either empirical or theoretical, of the validity of these models in the very low dose range. On the other hand, the one-hit model does represent an attempt to construct a rough causal model for low doses and hence many researchers studying radiation as a carcinogenic process favor the one-hit model of the process while others theorize “that it takes some minimum number of molecules, not a single molecule, to produce an effect” (Schneiderman (1974) p. 72)). The theoretical analysis by Arley (1961) and the analysis of virus data by Bryan (1961) are just two of many studies supporting the one-hit and the multiple-hit models in the case of carcinogenesis. Yet in their present form, extrapolation based on these models leads to extremely conservative results and hence they are not widely used. The ~ompa~mental safety model herein proposed and also previously investigated by Matis, Cardenas and Kodell (1974) generalizes the theoretical appealing one-hit model to a multiple-hit model within the broader framework of compartmental analysis. Compartmental analysis is a technique which has broad acceptance in the biomedical sciences in general and in physiological modelling in particular (see e.g., Atkins (1969), Jscquez (1972), Rescigno and Segre (1966), or Sheppard (1962)). A compartmental model typically does not purport to be a detailed causal model of a phenomenon but instead represents a rough and readily eon~ptualized framework su~ciently rich in parameter to provide an adequate approximation to many physiological systems. Moreover the recent addition of stochastic variation to the compartmental model (Matis and Hartley (1971)) provides a rationale for making valid inferences about a complex phenomenon from the relatively simple compartmental model
USE OF A STOCHASTIC MAMMILLARY
COMPARTMENTAL MODEL
469
inasmuch as the need for the intricate detail of the phenomenon is supplanted by probabilistic variation. The theoretical justification of the one-hit and multi-hit models combined with the general acceptance and versatility of compartmental models make a strong case for the present safety model. Although most applications of the cited safety models deal directly with carcinogenic safety testing, the basic principles apply to a wide variety of other environmental and ecological applications. It should also be pointed out that there exists another class of safety models based on a time to response (tumor) relationship (see e.g., Armitage and Doll (1961), Peto and Lee (1973), and Hartley and Sielken (1975)) which is not considered in this paper. 2. ilfodel I-An Upen M~~nrn~ll~T~ Compartment 3fodel. We assume the mammillary compartmental model illustrated in Figure 1 and designated as Model I. Let Xi(t), for i = 1, , . . , k, denote the number of particles in each of
h. a
1
Figure 1.
K2=
Ia x,(t)i 1”:
r5 a
.
l
0
x,(t)
A simple mammiIl&y comp~rtmen~l compartments
‘c20
i
model w&h k peripheral
the peripheral compartments at time t and Xk+l(t) the number remaining in the central compartment. It is also assumed that the central compartment is pulse lnhelled with a known number, say N, of independent particles at time 0, \vhereupon Xk+i(O) = M and X,(O) = 0, i = 1, , . ., k. The parameters a, ~1, and 7~sdetermine the fixed transition intensity c~~cients according to the following definitions (for small At): Prob. Prob. Prob.
{a, given particle in camp. k + 1 at time t transfers into camp. i in the interval i = 1, . . ., k. (t, t+At)) = a&+0(&), {a given particle in camp. k+ 1 at time t leaves the system in the interval (t, z+At)f = x,aA.1f O(At), {a given particle in camp. i at time t leaves the system in the intervai (t, t+ At)] = K2aAt+ OfAt), i = 1, . . ., k.
These assumptions on t,he transfer p~obab~l~t~es are analogous to the exact transfer proportions specified in t.he more restrictive deterministio model. ,The model is termed “open” in the peripheral compartments for ~2 > 0.
470
J. H. MATIS, R. L. KODELL AND JI. CARDESAS
In a physiological safety model, the central compartment may be envisioned to be t,he blood plasma and the k peripheral compartments are either cells or groups of cells. This conceptualization is suggested by the Jacquez et al. (1960) model of the distribution of drugs in the body. A reaction (e.g. carcinogenesis) may be assumed to occur in the model if at any time any of the k peripheral compartments attains a specified threshold size, say of s particles. The notion of both fixed and stochastic thresholds has been widely used in carcinogenesis, and the above reaction would be characterized as an s-hit, single stage mechanism by Neyman and Scott (1966). Finally, the openness property allows repair within the peripheral compartments. As an example of an environmental application, a cotton field represented by the central compartment and the k adjacent fields by the peripheral compartments. An adjacent field becomes very attractive to and hence also infested with boll weevils if ever a threshold number of s such insects chance to be present at any time. Whatever the application, Model I undoubtedly represents a gross oversimplification of the phenomenon and many generalizations are possible, some of which will be discussed subsequently. Nevertheless Model I as presently constituted has the great advantage of parameter parsimony and yet at the same time provides much richer models for extrapolation than those presently available. It was shown in Matis, Cardenas and Kodell (1974) that the multivariate probability distribution of the random vector Xi(t), X,(t), . . . , Xk+l(t), X,+,(t), where k + 2 denotes the exterior compartment, is a (k + 2)-variate multinomial with parameters N, pi(t), . . . , pi(t), pz(t), pa(t) where pi(t) = (k + ~1 - Ks)-l{exp ( - Kzat) - exp ( - a(k + p2(t) =
exp (- a(k+ K$),
~30)
1 --h-W
=
KI)~)),
and
-pzV).
Letting Xi denote the maximum which camp. i ever attains over time, and Y’ the maximum of Xi, Xk , . . . , XL, it follows that Prob [Y’ 2 s] denotes the probability of a reaction. The analytical expression for Prob [Y’ 2 s] is intractable but it is shown that a lower bound is Prob [Y’ 2 s] 2 1 - {Prob [Xi < s]>~ 2 1-
8-l
c
{ f-0
jy 0 i
where p = (K~)-~c-~‘(~-~) with c = (k + K~)/KZ . Although the above is useful in its own right, a safety model requires an upper bound on Prob [Y’ 2 s] for use in subsequent extrapolation. We propose the procedure below for generating this upper bound.
USE OF A STOCH.4STIC
MAMMILLART
COMPARThfENTrlL
MODEL
471
3. Model II-A Closed Mammillury Compartmental Model. The closed mammiilary comp~tmental model, where ~2 = 0 and consequently where each peripheral compartment is a ‘%rap” from whence particles cannot escape is designated Model II. Lettin g Tg(t) denote the number of particles in compartment i, it follows that Cfz,s T,(t) = N. In the safety model context, a harmful reaction would be assumed to occur if any peripheral compartment reaches a threshold of s cumulative part,icles. Model II, in addition to providing an alternative safety model, also yields an upper bound for Model I as follows. Let T; denote lim T,(t) over time. Since the particles are independent and the eompart~nents mutually exclusive, it follows that the multivariate probability distribution of !l’; , Ti, . . . , l’L+z is also a (k + 2)-variate multinomial with parameters N, pr, ~1, . . . , pi, 0, p3; where pl = (k+ &-I and 13s = ~i(k+ ~$1. Let S’ be the maximum of T; , Th , . . . , Ti, i.e., the maximum peripheral compartment. One may then conclude from t,he Bonferroni inequality (Feller (1968)) that Prob[S’
fe,
h-2,8).
472
J. H. MATIS,
R. L. KODELL
AXD
M. CARDENAS
The entries in Table I illustrate numerically the performance of the bounds, fa and f4, as approximations to the true probabilities, fi and for specified thresholds, 8. Values offi and fz were generated by computer simulations with a = 1, k = 10, rcl = 10, and ~2 = 1 or 0 for Models I and II, respectively. These empirical values offi and f2 are the result of 1000 realizations of the given models for each value of N. The values of f3 and f4 were determined from expressions (1) and (2). For a physiological safety model, the application mentioned above becomes inverted. That is, for a specified, typically very low risk of p*, the low-dose
fz,
TABLE I Numerical Illustration of Probability Bounds N
40
60
80
100
S
f2
fs
f4
fl
5 6 7 8
0.2392 0.0648 0.0139 0.0025
0.341 0.087 0.019 0.003
0.413 0.124 0.030 0.006
0.4803 0.1388 0.0339 0.0071
6 7 8 9
0.3527 0.1297 0.0383 0.0097
0.496 0.204 0.057 0.012
0.612 0.280 0.094 0.027
0.7872 0.2970 0.0979 0.0285
8 10
0.1908 0.0688 0.0214
0.282 0.099 0.034
0.416 0.163 0.059
0.4669 0.1840 0.0654
9 10 11
0.2433 0.1008 0.0363
0.356 0.155 0.058
0.506 0.254 0.111
0.6309 0.2819 0.1147
9
“safe” level of Model I, say nz , or of Model II, say nf , are required, but they are available only through extensive simulation. In their stead, a “conservative” dose nl,could be derived from the upper bound function, f3, and the lower bound function, f4, could be used to monitor the sharpness of the upper bound, since f3(nz)-f4(ng) represents the maximum possible deviation of f3 from the true probability, be it fi or f2. If, on the other hand, the application required a “liberal” dose, the roles of f3 and f4 would be reversed, and a liberal dose, nz, could be derived. Since, clearly, nj s ni 2 n: 2 nz, the difference nz - ng would also provide a convenient monitoring device, whatever the application. The interrelationships of fi, f2, f3 and f4 as functions of iV are represented graphically in Figure 2. Input parameters were those used to generate Table I,
USE OF A STOCHASTIC MAMMILLARY COMPARTMENTAL MODEL
473
the threshold level being s = 8. There appears to be good convergence in the low-dose extrapolation region of the curves, as will be borne out by the following application. Suppose that in an environmental problem, values of s = 8, k = 10, ~1 = 10 are known parameters relating to a certain chemical agent for which a low-dose “safe” level is desired. Assume further that a risk level of 1 in 10,000 haa been ascertained as acceptable. Thus, regardless of which of the two compartmental
DOSE N
Figure 2. The probability of exceeding the threshold under Model I (fi) and under Model II (fi) with upper and lower bounds (A andfa) as a function of dose size N. The parameters for this example were 8 = 8, k = 10, a = 1, KI = 10 and Kz = i or 0
models applies, to determine the value of nl corresponding to a risk level of p* = l/10,000, one solves f&s)
= 10{ 1- ,-,(y)(-L)” i
(LY)“-‘1 = 0.0001
for N = ni . Repeated evaluation offs(N) on a computer for a range of N-values yields f3(23) = 0.00009713 and f3(24) = 0.00013928. A conservative choice of dose for the given risk !evel would then be ng = 23. A check on the precision of
474
J. H. MATIS,
R. L. KODELL
&SD
JI. CARDESAS
the probability bound leading to this choice of close is made by computing f4(23) = 0.00003042. An appreciation of the excellent capability of the extrapolat8ion procedure can be gained if the acceptable risk level is decreased drastically. Consider now an acceptable risk level of only 1 in 100,000,000, a risk level which is common in the safety testing of potential carcinogens. To determine an acceptable dose, ng ,one solves fs(N)=l
$ I-
,_,( ; X&.~(~)“-“)
=0.00000001
for N =nz. As before, Orialand error yieldsfs(9) =3.359x 10Tgandfs(10) = 1.605 x lo-*. Again, a conservative choice for an acceptable dose is nz = 9. The corresponding lower bound for a close of AY= 9 and a threshold level of 8 = 8 is f4(9) = 9.581 x 10-l’, thus indicating the precision of the upper bound on risk. 5. On the Implementation
of the safety Node1
The previous trial and error procedure involved in finding the safe dose, nz, by solving for N in p* =fe(N) can be simplified somewhat by using an initial estimate of 1~: based on an approximation to fa. Let a new constant E.be defined as (a) Initial
approximation
of nz.
#I = l-
p*/k =
s-1 iv 1 . ti=o0 a
pf(1 -p)N-“.
Inasmuch as p is small, one can use the Poisson approximation to the binomial to obtain R N CfZi exp (-Np)(N p )i/‘l z ., w h’lc h CZIAbe rewritten as the cumulative ~2 distribution (see e.g. Beyer (1963)) as A=
s
O”exp
( - x2/w2/w-1
dX2
2r(s)
~NP
Hence 2Np is the upper 1 percentage point of a ~2 distribution with parameter 29, usually denoted x&r+ and therefore nz 2: xzs, 142~. To illustrate the approximation, assume again that s = 8, k - 10,~1 = 10 and p* = 0.0001, from whence p = 0.05, 1 - i, = 0.00001 and n3 =
(x:3,
o.oooolW~.
Using Table II by Lentner (1975) for the extreme lower percentage points of the ~2 distribution, one has ns = 1.9931/0.1 N 20. With this initial estimate, the procedure would rapidly converge on the safe dose nz = 23. In the same
_
_
0.05233104 0.13020332 0.25604775 0.43368150 0.(16383980 0.94544407
0.~000000 0.00002828 0.00168722 0.01401817
0.~~~1
0.08315227 0.19195597 0.35804763 0.58303141 0.86606411 1.20426967
0.00000000 0.00008944 0.00303589 0.02495551
0.00000000~ -
0.13232794 0.28359924 0.50198774 0.78641547 1.13355429 1.53921731
0.00000002 0.00028286 0.00783741 0.04446447
0.000~00~
_F
~I__-
--.-
._
_
2-‘qI-(v/2)1-1
0.22110412 0.42033933 0.70644114 1.06496937 1.49021874 1.97637282
0.00000020 0.00089450 0.01690433 0.07934626
0.0~~01
t Courtesy of Marvin Lentnor, Dept. of Exp. Stat., New Me?tioo &da Univ.
10 12 14 I6 18 20
2 4 G 8
28
II
Percentage Points of ~2 Distribution* F[x~/v] =
TABLE
0.25535777 0.54516954 0.93959477 1.42018445 1.99310094 2.F3012172 3.32867533
0.14195478 0.33812600 O.F2GO(i29G 0.99903807 1.45069332 1.97113791 2.55363768
0.00001 ___
2/2)aY2-1
0.0~002000 0.00895763 0.07907433
-~
(-
0.000002~~0 0.00282970 0.03G50857
0.000001
j”o* exp
ckc
0.40359379 0.88892030 1.42749511 2.06079786 2.77393852 3.55516529 4.39516272
0.00020001 0.02841848 0.17235211
O*OOOl
476
J. H. MATIS,
R. L. KODELL
AND
Bf. CARDEXAS
problem with p* = 0.00000001, it follows that ns = Cxz, o.ooooOOoOl)/O.l = 0.5330/0.1 = 6, as compared to ~2.:= 9. It can be shown that an approximation to 1~: is also -C/C- 1 . 124= x&+JZ~ with 1 = (I- p*) ’ and p = (tc#1 c (b) linknown parameters k, ~1, ~2, 8. In many environmental applications the parameters have a well-defined structural interpretation and they are either known or measurable by suitable ex~rimentation. In physiological applications, however, the parameters may not correspond as closely to recognizable physical characteristics and hence it is less likely that the parameters can be specified initially. As with other safety models, the case of unknown parameters necessitates testing the response at multiple dose levels, say where ~lt exceeds the number of u~no~ parameters. Xl, N2, . . . , N, Equation (I) or (2) may then be fitted to the dose-quanta1 response data and the parameters thereby estimated. Although the estimation procedures raise some auxiliary statistical questions, the same general extrapolation principles hold. The analytical study and simulation of the statistical procedures are a topic of current research. (c) Generalizations. It is not difficult to postulate a multitude of realistic generalizations of this basic mammillary model. One immediate generalization would arise by allowing unequal transfer coefficients both into and out of the ~ompa~ments, say olg and 66, i = 1, 2, . . . , k, respectively, as considered previously by Matis, Cardenas, Kodell(1974). Another natural extension would be to incorporate a stochastic threshold where s would have a probability distribution and one could also plausibly argue that the transition intensity coefficients must almost surely be time dependent and/or that N is a random variable with a known probability distribution (see e.g. Cardenas and Matis f 1974)). Although involved expressions could be written for each of these generalizations, it is likely that the present model is very robust against such alternatives. However as data in the low-risk dose range becomes available, the validity of the basic model will be tested and its structure refined to fit particular applications. 6. Discussion. The incorporation of the multiple-hit reaction mechanism within a compartmental structure suggests a broad class of safety models of which the present safety model using a particular mammillary system is only one special case. Besides their conceptual appeal, those models offer a family of curves for extra~lation from the high-risk experimental doses to the low-risk use doses. It may also be possible to extrapolate from “mouse-to-man” with these models by correlating a set of species specific parameter estimates with
USE OF A STOCHASTIC MAMMILLARY COMPARTMENTAL MODEL
477
some well-defined physiological characteristics, thereby extrapolating from a hierarchy of animal species to man. The referees of our previous paper, R. Elashoff, A. Rescigno and C. W. Sheppard, made several suggestions for its extension, some of which are incorporated in this manuscript. We acknowledge and appreciate their comments on the previous and the present paper. The investigation was supported in part by a Public Health Service Research Career Development Award 1 K04 BM00033-01 from the National Institute of General Medical Sciences to the first author.
LITERATURE Arley, N. 1961. “Theoretical Analysis of Carcinogenesis.” Proo. Fourth Berkdey Symp. Math. Stat. Prob., 4, 1-18. Arm&age, P. and R. Doll. 1961. “Stochastic Models for Careinogenesis.” Proc. Fourth Berk~y Symp. Math. Stat. Preb., 4, 19-38. Atkins, C. L. 1961. ~~~~t~orn~~n~ent Models in Biological Sy8teme. London: Methuen. Beyer, W. H. (ed.) 1968. Handbook of Table8 for Probability and Stati8tie8, 2nd Ed. Cleveland: Chemical Rubber Co. Bryan, W. R. 1961. “Virus Carcinogenesis.” Proc. Fourth Berkeley Symp. M&h. Stat. Prob., 4, 123-152. Cardenas, M. and J. H. Matis. 1974. “On the Stoch&ic Theory of Compartments: Solution for ~-Compa~ment Systems with IrreversibleTime-De~ndent TransitionProbabilities.” Bull. Math. Biol., 36, 489-504. FDA Advisory Committee on Protocols for Safety Evaluation. 1971. “Panel on Carcinogenesis Report on Cancer Testing in the Safety Evaluation of Food Additives and Pesticides,” Toxic. Appl. Pharmac., 20, 419-438. Feller, W. 1968. An Intmduction to Probability Th,ewy and It8 AppEieatiaow, Vol. 1, 3rd Ed. New York: Wiley. Hartley, H. 0. and R. L. Sielken. 1975. “A Parametric Model for ‘Safety’ Testing of Carcinogenic Agents.” FDA Technical Report No. 1, Institute of Statistics, Texas A 6 M University. Jacques, J. A. 1972. Compartmental Andysia in Biology and Medicine. Amsterdam:
Elsevier .
R. Bellman and R. Kalaba. 1960. “The Distribution of a Drug in the Body.” B&E. Xath. Biop~ys., 22, 309-322. Lentner, Y. 1975. Percentage Points of x2 Distribution. Personal Communication. JIant,el, N. and W. R. Bryan. 1961. “Safety Testing of CarcinogenicAgents.” J. Natn. CancerInst., 27, 321-346. Matis, J. H. and H. 0. Hartley. 1971. “Stochastic Compartmental Analysis : Model and Least, Squares Estimation from Time Series Data.” Biometric& 27, 27-102. &f. Cardenasa.ndR. L. Kodell. 1974. “On the Probability of Reaching a Threshold -7 in a Stochastic ~~a~ill~ System. ” Bull. Hath. Bid., 36, 445-454. Neyman, J. and E. Scot,t, 1966. “Statistical Aspects of the Problem of Carcinogenesis.” Proc. Fijth Berkeley Symp. Math. Stat. Prob., 4, 745-776. Feto, R. and F. Lee. 1973. “Weibull Distributions for Continuous Carcinogenesis Experiments.” Biometrics. 29, 457-476. Rescigno, A. and C. Segre. 1966 l?r~cgand Tracer Kinetics. Waltham: Blaisdell.
-,
478
J. H. MATIS, R. L. KODELL AND M. CAFtDENAS
Schneiderman, M. 1974. “Symposium Acud. Science. Sheppard, C. W. 1962. Baeic Primiples
on Statistics
and the Environment.”
J. Wash.
of the Tracer Uethod. New York: Wiley. RECEIVED 7-28-75 REVISED l-8-76