Power laws in compartmental analysis— II. numerical evaluation of semi-Markov models

Power Laws in Compartmental AnalysisII. Numerical Evaluation of Semi-Markov

Models*

ALLAN H. MARCUS AND

ARTHUR BECKER University of Maryland-Baltimore

County, Catonsville, Mayland

Received 30 September 1976

ABSTRACT A semi-Markov model for tracer kinetics is formulated for certain two- and three-compartment systems for which the residence time in compartment 2 may have a stable or exponentially modified stable distribution. The method requires numerical inversion of complex Fourier transforms. The model is used to predict the proportion of a single dose of lead retained by blood, bone, and soft tissue, where the lead residence time in bone has a distribution of power-law type. The power-law models reproduced the sometimes observed increase of blood lead level due to slow release of lead from bone, whereas the usual linear compartment model with matched parameters did not show this effect.

INTRODUCTION In Part I [4] we considered a unified random mechanism, based on the expected number density of residence times of tracer molecules attached to a random number of sites within a compartment, that produces mixed power-law exponential distributions of residence time, including not only Gamma distributions other than exponential, but also exponentially modified stable laws. We will now develop some techniques that can be used to model organisms consisting of several such interconnected compartments. The functions to be computed include the amount of tracer excreted in urine, feces, etc. at any time t, as well as the proportion of initially injected or ingested material residing in any other specified compartment at time t. *This work was supported by the National Institute of Environmental Health Sciences and the Environmental Protection Agency under Grant l-ROl-ES01236-I. The University of Maryland Computer Science Center provided additional computer time. MA THEMA TICAL BIOSCIENCES

0 Elsevier North-Holland, Inc., 1977

35,27-45

(1977)

27

28

ALLAN

H. MARCUS

AND ARTHUR

BECKER

The complement of the total quantity excreted is of course the remaining whole body burden (retention), and the proportion in any specified compartment is the body burden of that compartment. Formally, we will regard the specific activity in excreted material as an empirical estimate of the first passage time from the initial compartment to the external environment. The relative body burden in any compartment will be regarded as an estimate of the probability that a molecule is found there at time t after insertion. The relative activity in the excreted material will be described by a first-passage-time (FPT) probability density Jo(t). The relative body burden at time t in compartment j, per molecule of tracer inserted into compartment i at time t = 0, will be described by the transition probability b,(t). These quantities are easily calculated using transform methods in the case that the tracer molecules move independently of each other according to a semi-Markov or Markov renewal process whose theory we shall now develop in detail. SEMI-MARKOV

MODEL

Consider a single molecule of tracer. Molecules of tracer may be bound to other molecules-e.g., carbon monoxide may be bound to hemoglobin in the blood, but also to other hemeproteins in the liver or kidney. Then, the transition to other compartments may be governed by the size and other physical-chemical properties of the carrier-i.e., both the residence time and the subsequent transition depend to some extent on the carrier. Define Pii = Probability that a molecule travel to compartmentj

i will next

in compartment o< Pg< 1, (

c

Pq=l),

j-0

G, (t) = Prob{ molecule resides in i for time d t, given that molecule starts in i and travels next to j}, F, (t) = Prob{ molecule

delay

J,(t)=

that starts in i reachesj

for the first time after


(2)

fF,W>

A0 (t) = excretion

(3)

or clearance

b,(t) = Prob{ molecule

delay = t}.

(1)

function

at time t,

(4)

that starts in i is found in j after a (5)

POWER LAWS

29

Note that

(6)

biO(r)=/o'LO(U)du~

Note that&(t) is a first-passage-time (FPT) probability density function. If transitions are independent of each other, governed by a stationary Markov chain with transition probability matrix [ = (P,), and if the residence time in compartment i (knowing that when we exit from i we next go to j) has distribution G,(t) and is independent of the duration or direction of all previous transitions, then the process is a semi-Markov process. It is not hard to show [3] that, if we define the core matrix

~(t)=(Z=,G&))=~CW(t)

(7)

and the FPT matrix

(8) then the Laplace-Stietjes C*(s)=

transforms

are

( (p”dq(t) = c*(s)]-‘{ [z-

(9)

Pg~me-“dGti(r))

0

F*(s)=

0

c*(s)[ I-

C”(s)

]pJZ}

-I.

(‘0)

Here Z is the identity and @ is the Schur product of two square matrices of equal size, i.e., for any such A, B, A 8 B = (a&).

Moments

are easy to calculate.

If the following

mii=

J0mrdGg(~),

so=

s

“t2dG&),

0

‘k,=

I mt2dFu (t), 0

are finite:

ALLAN H. MARCUS AND ARTHUR BECKER

30 then N by=

*,=

z k=O

Pikmak

$

Piksik+

k=O

+

2 k=O k#j

5

Pikpkjf

(11)

Pik[9kj+2mik(1,,].

(12)

k-0 k#j

The body burden in compartment j at time t for a unit impulsive insertion of tracer into compartment i at time 0 is denoted bti(t). It is clear that b,,(t) can be calculated from a renewal-type integral equation likefii(t), so we find [3] the Laplace transform of the matrix

is given by B*(s)=/~Pb,(t)dt=

[I- c*(s)] -%v*(s),

(13)

0

where W*(S) is a diagonal

matrix with

(14) See [1,2] for an extended discussion of semi-Markov models. These functions can be calculated using inverse Fourier transforms. example, for a real variable z we can write

For

f:O( - iz) = Ref( z) + i Imf( z),

(15)

b$(-iz)=Reb(z)+iImb(z),

(16)

where i = (- 1)‘/2, and Ref, Im f, Re b, Imb are real functions. Note that Re,f and Re b are even, and Imf and Im b are odd functions of z, so that we need only evaluate the real integrals

for sufficiently

fi,,(t)=

$iAIRe_f(z)coszr+Imf(z)sinzt]dz,

(17)

b,(t)=

~~BIReb(z)coszt+Imb(z)sinzr]~z,

(18)

large A, B. The remainder

of this paper

will consist

of

31

POWER LAWS

exhibits of the functions Ref, Imf, Reb, Imb and the appropriate values of A,B. For reasons of computational simplicity we choose the following cases in which to exhibit explicit integrands. EXPONENTIAL

DISTRIBUTION

g,(t) = kuemkut,

t >O,

=o,

t GO,

g;(s) =

&.

(19)

rJ

POWER-LA go(t)

W DISTRIBUTION

OF STABLE

is not usually available, g$(-iz)=exp{

TYPE

but g;(s)=

e-y@,

so

-yi/jzI”[cos(ncu/2)-isgn(z)sin(~a/2)]}

or g$( - iz)=exp[

- y,ijzI”cos(acr/2)]

{ cos[ yJzlasin(sa/2)]

(20)

+isgn(z)sin[y,i(zl”sin(~cu/2)]}, where k,>O, THE

SPECIAL

yo>O, O
MODEL

As a result of straightforward but tedious calculations case depicted in Fig. 1 (compartment 0 leads)

we obtain

for the

so that in particular

and l-c*c*12 s[ 1 - c;~c&]B*(s) = i

21

0

0

47

I-cf,-c$

cz:c:o

c&(1 - Go-

42)

cr2(l -4) l-c,:

, I

(23)

32

ALLAN

H. MARCUS

SPECIAL TWO-COMPARTMENT

FIG. 1.

Special

AND

ARTHUR

BECKER

MODEL

two-compartment

model.

so that in particular 1 - CQ - it) - c&( - iz) bT,(-

iz)= i

’

z[l-cfi(-iz)c;,(-iz)]

b&( - iz) = ic$( - iz) This special two-compartment strictly internal and excretion

1 - cF2( - iz) z[1-cc:2(-iz)c:*(-iz)]

(24) (25)

model is one in which compartment 2 is can occur only from compartment 1. Thus,

p2l.l= 0,

pfl= 1,

P,*= 1 - P,,. The special cases are easily developed. First case: Al2 residence times are exponential. Define

A, (z> =

z4+ z’( k:z+ kt, + P&,&2,) +P&&I,

A2 =P&&,

(k12+ k21).

(26) (27)

Then f$( - iz) = Ref,, + i Imf,,,

(28)

Refi,(z)=Q,(z)[k,,A,(z)-2*A21, Imfio(z)=Ql(z)z[A~(z)+k~~A2],

(29) (30)

33

POWER LAWS

where (31) (32) (33) (34) where I

Q2(Z)= (k$)+z2)[A:(z)+z54;] h(Z)

=

(35)

’

z4(plob3+p12~12) +

z’[

klok,,(p,okl2+pl,k,,)+

+

k,&,,k~l

k,Zl (p,&,o+p12k12)]

(36)

(plokl2+pl2klo)~

~4(z)=~4+~2(k:0p12+k:~p1~+kf2) +

kt,

(p&:2

(37)

+pl2%,

b~2=Reb,2+iImb12,

(38)

Reb12=Q~(z)[Al(z>(k12k21-z2)-A2z2(k~2+k2~)],

(39)

Im~12=zQ3(z)[A1(z)(k12+k21)+A2k12k2~]l

(40)

where (41) Second case: Residence time in 2 has a stable distribution. Define L,(z)=exp(

- ,,2~zl~cos~)coS(~12/zJaSin~),

(42)

L,(z)=exp(

- y,,[z("cos~)sin(y,,(zlysin~),

(43)

H, (z) = k&[ 1 -~ltb

]

+

z2+p12k121zl~2(z)~

(45)

~2(~)=p12k12[~~l(~)+~gn(~)~2(~)k12]~ ff3

(z)

=

k12sdz)~2(z)

ff4

(z)

=

kl,[

(20(z>=

1-

W)]

(44)

+

1 H: (z) + H,2 (z) ’

z[

1-

lzlLz(z)~

Ll

(z)],

(46) (47) (48)

34

ALLAN H. MARCUS AND ARTHUR BECKER

~mk,o(k:,+z*)Q&) Q~(z)= kfo+z2 ’ Q,(Z)~

(49)

e, k$,+z2

(50)

’

Q6 (z) = Q,(z)P,~/z>

(51)

A5(z)=z2(~,0k10+~,2k12)+k,0k,2(~,0k,2f~,2k,0),

(52)

AC.(Z) = z(&:o+W%)+

z3.

(53)

Refio=Q4(z)[k10Hl(z)-zH2(z)],

(54)

ImS,,=

(55)

Then

Q4(z)[zH1

(z)+k~H~(z)],

Re~~~=Qs(z)[H,(~)AS(z)+HZ(z)A6(z)],

(56)

Imbll=

Q5 (z)[ H, (Z)&(Z)+

(57)

Reblz=

Q~(z)[H,(~)W~(~)+H~(Z)H~(Z)],

K (z)~,(z)],

(58)

Imb12= Q&)[Wz)&(z)+~2(~)%(z)].

THE SPECIAL

THREE-COMPARTMENT

(59)

MODEL

The case depicted in Fig. 2 is one in which compartment only, but excretion from either 1 or 3 is possible. Thus p,,+

2 is internal

p,z+ p13= 1,

(60)

Pzl+P23=

1,

(61)

P30= 1,

(62)

Poe= 1.

(63)

p,,+p32+

In view of the complexity of the special two-compartment report only the transforms of greatest interest:

model we will

A = (1 - c;c& - c&c$ - 4~3: - c$c&c;, - c:~cz;c;,),

bT2=

(l-c:,-c~3~(c:,c:i+c:,) 1.

(67)

Z

These results are readily deduced

(64)

by a standard

flowgraph

analysis.

Fig. 2

35

POWER LAWS

SPECIAL THREE-COMPARTMENT

FIG. 2.

STRATEGY

FOR NUMERICAL

MODEL FOR LEAD

Special three-compartment

model.

ANALYSIS

The inverse Fourier transforms were calculated by numerical integration of the Fourier transforms using Eqs. (17) and (18). We found that the positive and negative contributions during one period of the harmonic component nearly canceled. Thus zeros of Ref(z)

coszt + Imf(z) 7r

sinzt = Int(z)

(68)

= Int(z)

(69)

and Reb(z)coszf+Imb(z)sinzt ?r

were calculated by binomial search. The periodicity was nearly 27r/ t after the first 20 or 30 zeros. The positive and negative contribution was integrated between each pair of adjacent zeros, then summed over pairs of zeros. The upper limits A,B were increased until the relative increment in the integral was less than some specified value. The programs are divided into three main phases: (1) Integration (2) Integration for zeros. (3) Integration amplitudes

of the first zero. through subsequent

zeros in pairs, with binomial

search

from zero to zero, with an assumed spacing n/t once the are changing only very slowly with increasing z.

The user may control

the length of phase (2) before entering

phase (3).

36

ALLAN H. MARCUS AND ARTHUR BECKER

Phase (1). The lower limit of the integral number z. from 1O-5/~ to 10e3/t, because power-law residence time in the system. z until Int(z) is negative. The first zero of search; call it zi. The first increment (usually is thus calculated by Simpson’s rule as

is taken as some small positive b;(O) may be + 00 if there is a is increased by steps of 0.01 /t Int(z) is located by binomial the largest) in computing b;(t)

LnInt(z) dz, s =“-I accumulated

until In Int(z)dz<(fraction) I G-1

The fraction

is usually

‘“-‘Int(z)dz. I 10

(71)

small, 10e4 to 10d5. If we have in addition lz,-zz,_,--7rP/fl<(fraction)r/t

to (71) (72)

for some other small fraction, we pass to phase (3), Phase (3). We stop searching for zeros and take

z,=z,-,+(7T/f), again stopping

(73)

if for some small fraction (74)

We did not calculate

b,,(t) directly by this method, lim b,,(t) t-cc

Thus 6To(- it) must necessarily origin and cannot be accurately the retention probability

since

= 1.

have a singularity of order inverted for large t. Instead,

R (t) = 1- &,(4

(75) - 1/iz at the we calculated

(76)

since R*(iz)=

/

meiZrR (t)dt

= -ol,iz-b~o(-iz)

r:

.f;cO(-iz)-1 iz

is easily inverted

numerically.

(77) (78)

POWER

31

LAWS

The three-compartment model was calculated numerically using standard complex integration subroutines from the Univac 1108 Math-Pack. These were verified by a numerical integration of the two-compartment model functions as real integrals for the special cases of the preceding section. Compartment 3 was exorcised by setting Pi, = P3j = 0 for all iJ# 3. In fact, this necessary numerical validation motivated the tedious algebra required to calculate Eqs. (26)-(59). We believe that these explicit formulae may be of independent interest and use to readers who do not have access to computers with readily available complex-function subroutines. The accuracy of these numerical results is not extremely great-we cannot really assert that the relative errors of the probabilities are less than about 0.1%. However, the results in Tables l-3 show that the calculated sums blo+ b,i + 6,, + b,, are satisfactorily close to 1.0 in value.

NUMERICAL

ANALYSIS

OF A THREE-COMPARTMENT

MODEL

FOR LEAD

Rabinowitz, Wetherill, and Kopple [6] studied the distribution and transport of a lead in a mature human male using the stable lead isotope Pbzm as a label. We will follow the basic themes of their analysis. The compartments are the external environment (i =O), blood (i= l), bone (i=2), and soft tissues (i=3). For transitions out of blood and soft tissue, we assume the usual linear

Probabilities

t Cm4

b,,(t)

0 1 2 4 6 12 24 36 48 60.1* 72 96 168 360

0 .28694 .45301 .6288? .71924 .83026 .87060 .87873 .88484 .89272 .89605 .90623 .93115 .96983

TABLE for linear (Exponential)

1 Compartment

b,,(t)

b,,(t)

1.OOOOO 0.63366 0.43334 0.24342 0.15859 0.58614(0.12500(0.43533( 0.25152(-2) 0.20446( 0.18994(-2) 0.16380(-2) 0.12965( 0.67895( -

0 .11270(.22864( .42621(.59082( .92403( .11526 .11628 .11206

*There was some difficulty

1) 1) 2)

b,,(r) 1) 1) 1) 1) 1)

2)

.10675

2) 3)

.10151 .91551(1) .67184(- 1) .29408( - 1)

in getting

an accurate

Model (i)

0 .68813(.91688(.85717(.63538( .19421(.16025(-2) .51413(-3) .44076(-3) .41905(-3) .39917(-3) .36077( .26651(-3) .11700(-3)

result

Sum

1) 1) 1) 1) 1)

3)

for t = 60.

1.OOOOO 1.00068 1.00090 1.00063 1.00045 1.00070 0.99996 0.99988 0.99986 1.00193 0.99986 0.99980 0.99990 1.00003

38

ALLAN

Probabilities

TABLE 2 for Exponentially Modified

f bo)

b,,(4

b,,(t)

0 1 2 4 6 12 24 36 48 60. la 72 96 168 360

0 .28696 .45302 .62874 .71886 .82840 .86373 .86602 .86694 .86913 .87284 .88393 .92187 .91489

1 .63286 .43227 .24243 .15756 .57213( .10915(.25486( .10473(-2) .11702(-2) .16737( .23261( .21263(.78190(-3)

‘There

H. MARCUS

was some difficulty

Stable-Law

1) 1) 2)

2) 2) 1)

0 .12182(.23419( .43288( .60173( .95630( .I2430 .13126 .13181 .12948 .12524 .11341 .75730( .24348(

-

1) 1) 1) 1) 1)

- 1) - 1)

an accurate

TABLE 3 for Stable-Law

ARTHUR

1 1.00080 1.00039 1.00012 1.00006 0.99994 1.00015 0.99992 0.99988 0.99997 l.OOW5 1.00012 1.00015 1.00015

result for r = 60.

Model (iii)

bd’)

b,,(O

b,,(r)

b,,(f)

0 .30321 .46508 .63162 .I2624 .83375 .86757 .86907a .86891

1 .63294 .43228 .24244 .I5756 .57255( - 1) .10619(-l) .24732( - 2) .98688( - 3)

0 .63179(-2) .19201(.40270( .57700( .90168( .12042 .12826 .13021

0 .68801( .91682( .85690( .63468( .18715(.12056(-2) .73996( .18612(-4)

60.1b 72 96 168 360

.86894 .86940 .87122 .87999 .90032

.60689( .50090( .70154(-3) .82679( .60030( -

.13043 .12995 .12792 .11907 .99799( - 1)

3) 3)

Sum

0 .68799( - 1) .91681(- 1) .85664( - 1) .63466( - 1) .18695(- 1) .12057(-2) .87260( -4) .84814(-4) .19132(-3) .30544( - 3) .45814(-3) .41989(-3) .13347(-3)

0 1 2 4 6 12 24 36 48

3) 3)

BECKER

Model (ii)

b,,(r)

b,,(t)

in getting

Probabilities

AND

1) 1) 1) 1)

Sum

1) 1) 1) 1) 1) 4)

.36602( - 4) .66517(-4) .10952(-4) .13226(-4) -

‘This suspicious value could not be calculated to greater accuracy able amount of computer time. bThere was some difficulty in getting an accurate result for t = 60.

1 1.01127 1.00824 1It0602 1.00497 0.99995 0.9998 1 0.99988 1.00900 l.OOWl 0.99992 0.99995 1.00002 -

within a reason-

39

POWER LAWS

compartment

model, i.e., a Markovian

model with

Pu= ku/ki, k,= 5

k,,

j-0

k,

g;(s) = kv+s ’ for i= 1,3. The parameters times

suggested

(81)

by [6] include

the mean

residence

may be deduced from the mass-balance

relation-

E{T,}=l/k,=27days, E { T,} = l/k,=30 The probabilities ships

required

from 1 to 0: from 1 to 2: from 1 to 3:

38 pg/day, 7 pg/day, 20 pg/day;

from 3 to 0: from 3 to 1:

12 pg/day, 20 pg/day;

from 2 to 1:

7 pg/day.

Thus (putting

days.

rates in units of months-‘) PI,= 38/65 = 0.584, p,2=

7/65=0.108,

k,,= 0.6489 moo’, k,,=O.l200mo-‘,

P ,3=20/65=0.308,

k,,=0.3422

P 3o= 12/32 = 0.375,

k30 = 0.375 mo- i,

P,, =20/32=0.625,

k3, =0.625 mo-‘,

p,,

=

mo-‘,

0,

where 1 mo= 30 days. In order to link the exponential and stable-law models to the bone data, we have had to include an intermediate type of distribution which was also derived in Part I, the exponentially modified stable law (EMSL). There are thus three models for the residence time in bone: (i) Exponential

model.

La(s)=

&

(82)

40

ALLAN H. MARCUS AND ARTHUR BECKER

whose

mean and variance

(ii) Exponentially

are

E { TmI= l/k,,,

(83)

Var{ T, } = l/k&.

(84)

modified stable-law model.

(85)

Y21>

k,,

0,

O
> 0,

Then

E { 7.2,>= w%,$3 21

a(l-a) Var{

r2,

> = ~2&2q ___. k:,

(iii) Stable-law

model (k,, = 0).

w

g2*I(s)=e-Y21So,

E{Tg,}=+co

for

p> a.

(89)

Now, Rabinowitz, Wetherill, and Kopple indicate only a very long residence time in bone. They suggest a mean residence time of about 10,000 days (28 years). They indicate nothing about the appropriateness or otherwise of the linear compartmental model (82) for bone. We will here study the implications of the hypothesis that, like the alkaline-earth metals calcium, strontium, barium, and radium [S], lead is a bone-seeker whose residence time in bone has either an EMSL or a stable distribution. To fit the numerical values, we will consider a somewhat shorter mean lifetime in the exponential case, 200 mo = 6000 days = 16.67 years. Thus for Model (i) we have I/k,,=2OOmo=E{T,,} 1/k& = 40,000 mo’ = Var { T2, } . We cannot simply equate moments to fit model (iii), the stable law, which has an infinite mean and variance. One possibility is the (approximate)

POWER

matching

41

LAWS

of modes: y;;

g21 (t)

=

(90)

g21 (Lde)~

where we easily find t mode =

(91)

~26.

For model (ii) (EMSL), (9 + y24rk:,)0~5- 3 t mode

=

YW,

.

(92)

Here we can set (Y= 0.5. Thus from (86) and (87) E { T,, } = 200 mo = 0.5~,,&‘.~.

We will here use the same exponential l/k,,

time scale

=200 mo,

so that y2, = 28.28 mo”,5, k,, = 0.005 mo- ‘, and from (92) t mode = 100 mo. Applying matching

the similar value fmode= 133 mo to (91), we find that the modevalue for model (iii) is also y2, = 28.28 mo”,5.

INTERPRETATION

OF THE LEAD

MODEL

We will describe elsewhere the detailed significance of these numerical results. In Figs. 3-6 we exhibit respectively the retention probability R (t) = 1 - bra(t), the blood level concentration b,,(t), the bone concentration !~,~(t), and the concentration !~,~(t) in soft tissues, on both semilog and log-log plots, for all three models given in the preceding section.

42

ALLAN

H. MARCUS

AND

ARTHUR

BECKER

WHOLE BODY RETENTION FUNCTION FOR MEAN LIFETIME 200 MONTHS 1

10-21 0

I

I

I

I

,

I

100

I

I

I

,

I

7

200

1

I

,

300

I

I

I

I

,

400

TIME. MONTHS WHOLE BODY RETENTION FUNCTION FUR MEAN LIFETIME 200 MONTHS 1-j

FIG.3. exponential, log-log

Retention function for lead models with residence time in bone given by (E) (EMS) exponentially modified stable, (S) stable distributions; semilog and

plots.

The retention probabilities exhibit no surprises, decreasing most rapidly for model (i) and least rapidly for model (iii). They are very similar to the retention curves calculated by Marshall et al. [5] for the alkaline-earth metals calcium, strontium, barium, and radium. The blood level concentrations b,,(t) shown in Fig. 4 are highly significant. Note that b,,(t) decreases monotonically for all times t, whereas with either stable or EMS laws, b,,(t) reaches a minimum after 30-70 mo, but climbs to much higher values after about 100-150 mo. We have thus reproduced mathematically the well-known inclination of blood lead levels to increase many years after the cessation of the last exposure. Only the power-law models did this; the linear Markov compartment model did not,

POWER LAWS

43 BLOODLEVELOF LEAD FOR MEAN LIFETIMEIN BONEOF 200 MONTHS

60.0 120.0 180.0 240.0 300.0 360.0 420.0

0.0

TIME.MONTHS BLOODLEVELOF LEADUSINGMEAN LIFETIMEIN BONEOF 200 MONTHS

10-Q

,

1

,

,

, , , ,,,,

10

,

,

, , , ,,,i

102

,

, ,

,,,flTIME.MONTHS 103

FIG. 4. Lead concentration in blood with residence times given by (E) exponential, (EMS) exponentially modified stable, (S) stable distributions; semilog and log-log plots.

and we conjecture that the methods used by Hearon [2] to establish the convexity or log-convexity of washout curves may be applicable here to establish this result in general. The bone concentrations b,,(t) in Figure 5 also show the expected behavior, peaking after 30-90 mo and decreasing slowly thereafter. Finally, the concentration 6&t) for the soft tissue compartment is shown in Figure 6. Notice the extremely pronounced increase in soft tissue concentration after 100 to 150 mo using either model (ii) or (iii), and not evident in the usual compartment model (i). There is also a modest short-term peak at 2-3 mo which is due to the release of a large quantity of the initial blood lead into soft tissue.

ALLAN H. MARCUS AND ARTHUR BECKER

44

BONE LEVEL OF LEAD FOR MEAN LIFETIME IN BONE OF 200 MONTHS

10‘3

I 60.0

c.0

1 120.0

I 180.0

I 240.0

I 300.0

I 360.0

, TIME.MONTHS 420.0

BONE LEVEL OF LEAD FOR MEAN LIFETIME IN BONE OF 200 MONTHS

lo-"

I

I

I , ,111,

1

I I 1,111

10

10-Z

I

t I I II~

TIME. I?ONTHS 103

FIG. 5. Lead concentration in bone with residence times given by (E) exponential, (EMS) exponentially modified stable, (S) stable distributions; semilog and log-log plots.

OTHER

SEMI-MARKOV

MODELS

A general semi-Markov formulation with many explicit results is developed by Weiner and Purdue [7]. An alternative formulation, also in the context of compartment models, is given by Wiggins [8]. Wiggins replaces the formulation based on residence time distribution by one in which the transition probabilities are taken as time-dependent; Howard [3] has shown that these formulations are equivalent. REFERENCES 1 E. Cinlar, Markov renewal theory, Ado. Appl. Probab. 1, 123-187 (1969). 2 R. B. Gun, C. S. Patlak, and J. Z. Hearon, The logarithmic convexity of the washout

45

POWER LAWS CONCENTRATION OF LEAD IN SOFT TISSUE FOR MEAN LIFETIME 200 MONTHS

0

300

200

100

400

TIME, MONTHS

10-S

I

,

, I I,

,,,

10

I

1

I I , ,111

102

1

I

I I I

,qT’NE.MoNTHS 103

FIG. 6. Lead concentration in soft tissue with residence times given by (EMS) exponentially modified stable, (S) stable distributions; semilog and log-log plots.

function in tracer kinetics, Math. Biosci. 4, l-6 (1969). R. A. Howard, Dynamic Probabilisfic Sysrems. Vol. II: Semi-Marka, and Decision Processes, Wiley, New York, 1971. A. H. Marcus, Power laws in compartmental analysis, Part I: A unified stochastic theory, Math. Biosci. 23, 337-350 (1975). J. H. Marshall et al., Alkaline earth metabolism in adult man, Health Phys. 24, 125-221 (1972). M. B. Rabinowitz, G. W. Wetherill, and J. D. Kopple, Lead Metabolism in the normal human: stable isotope studies, Science 182, 725-727 (1973). D. Weiner, and P. Purdue, A semi-Markov approach to stochastic compartmental models (1976), unpublished manuscript. A. D. Wiggins, A mathematical model relating the power law to the exponential law in biological turnover studies, Math. Eiosci. 10,191-200 (1971).