Compur.
Biol. Med.
Pergamon
Press 1976. Vol. 6, pp. 33-51.
Printed in Great
Britain.
NUMERICAL SIMULATION OF A COMPARTMENTAL SYSTEM FOR THE DISTRIBUTION OF LIPID SOLUBLE CHEMICALS IN MAMMALIAN TISSUES SHARON E. RODECAP Departments
and F. T. LINDSTROM
of Mathematics and Statistics, Oregon State University, Corvallis, Oregon 97331, U.S.A. (Received 23 October 1974)
Abstract-An analysis of the stability and accuracy of a particular numerical scheme used in the mathematical modeling of the flow limit pharmacodynamics of highly lipid soluble drugs in mammalian tissue systems is presented. The necessary definitions and theorems together with their proofs are included. A typical example of the flow limit pharmacodynamics of the highly lipid soluble drug dieldrin (HEOD) is also presented. Numerical dieldrin
Mathematical model
Flow limit pharmacodynamics
Lipid soluble drug
1. INTRODUCTION In recent years, with the development of high speed digital computers, the use of mathematical compartmental models to describe and predict biological phenomena has become more prevalent. One of the major uses has been to describe blood flow and drug distribution in mammalian tissue systems. A compartmental model is set up by partitioning the body into a finite number of compartments [ 11, each of which consists of tissues which are relatively homogeneous with respect to the drug being studied. To each compartment is assigned a variable denoting the amount of drug in that compartment. A system of differential equations is then developed giving the time rate of change of the amount of drug in each compartment in terms of the current amounts in each of the compartments, input from outside sources and excretion and destruction of the chemical. Correct simulation of these systems involves knowing a priori the compartment masses (volumes) and flow rates of the blood (in which the drug is assumed to be dissolved) between compartments. Using this basic method Lindstrom et al. [2] have formulated a preliminary model to describe the pharmacokinetics and the pharmacodynamics (drug distribution mechanics) of highly lipid (fat) soluble drugs in mammalian tissue systems. The following assumptions are made : (I ) The drug is transported via the blood stream completely dissolved in blood lipids (fats). (2) Drug transfer from the blood capillaries of the lipid fraction of each compartment is rapid (flow limit conditions operating). (3) The drug is converted (metabolised) only in the liver and is itself not easily excreted by the kidneys. 33
34
SHARON E. RODECAP and F. T. LINDSTROM
(4) Tissues can be partitioned into relatively homogeneous compartments within which there is uniform, instantaneous mixing (no intracompartment concentration gradients, difTusion currents, exist). (5) There exists a complete system of mass balance equations describing the drug distribution. (6) All blood lipid flow rates (QJ and compartmental lipid masses (M,) are continuous functions of time (0 I t I x). Note that these assumptions and the following mathematics are general enough to be useful for simulation of any highly lipid soluble drug in the body. Figure 1 is a diagram of the system with rectangles representing compartments, I/;
I
’
(Fecal excretion)
Fig. 1.
the amount of drug in compartment i (g), Mi the lipid mass (g), Qi the flow rate of blood lipid through the ith compartment (g/hr), 9 the input (source) function (g/hr). and r the liver drug conversion rate coefficient (g/hr). The remaining letters represent various biological constants. The lumen, although shown as a rectangle, is not a regular compartment; but it is simply a part of the liver compartment into which the drug enters with some slight excretion of the drug allowed. The term “mass balance equation” in assumption (5) simply means h:
the time rate of change of drug mass dnJdt i
iI
= [drug mass flow delivered via arteries from the general blood pool] + [drug sources] - [drug conversion, (sinks)] E) ( - [drug mass removed from compartment via venous blood]
Numerical simulation of a compartmental
system
Using this basic concept of mass balance on each compartment, equations is derived:
&L dt
0 .y + S?QH,f$ z
du,
z=~,z+Qs
-
8 du, dt = Q32
!+
Q,
2 -
Q,
35
the following system of
nil;-, (G.I.),
x
(aQHI + Q, + Qz + r)$
- (Q3 + R,,)$
(1)
I
Z
(Liver),
(Kidney),
(2) (3)
3
$f
=
Q4s
x
- Q3$
Q5$
- Q,$
du, dt - Qh 2
- Qh$
2
=
du, dt = Q7 $
x
(Brain),
(4)
(Depot fat),
(5)
(Muscle),
(6)
(Other tissues),
(7)
4
- Q, $
du, dr - (Q, + Q,,$
2
5 h 7
+ i ai .j=3
- Q8$
(Blood pool), 8
I
(8)
with the auxiliary equation: 7
Q8 = z
Qj (Conservation
of
flow).
j-l
(9)
Equations (l-8) can be written - = AD + s(t). dt
(10)
where OF-
I:_ 0
0
S=
.
9
0
and the matrix A appears in Fig. 2. Note that all the mass and flow parameters Mi and (2; respectively, are assumed to be real positive piecewise continuous functions of time t over CO,7J, T> 0.The special terms, 13, R,,,, 9, and r are positive real parameters with 0 I 8 I 1. The term 9. the source function, is assumed to be at least bounded and measurable [0, 13. For the most part R,,,is simply assumed equal to zero.
SHARON
E.
RODECAP
and F. T.
LINDS~OM
-Q1
-5 Ql
T -(Q3+Rexl 0
0
0
M3 0
0
0
0
0
0
0
0
0
Ql+Qz --
-Q4 M,
0
0
0
0
0
0
0
“b M6
0
0
0
0
-Q5 M5
0
0
0
0
M2
Fig. 2. Mammalian transport matrix for a male mammal.
2.
FORMAL
SOLUTION
Consider the first order system dU - = AD + S(r), 0 I t < cx, dt D(O) = 0, det(A) # 0,
with a;,i = aij(t), i= 1,2,. . . ,r; j = 1,2,. . . , Y continuous [O 5 t < Y&T > 0) and with all elements of S piecewise continuous on the same interval. This can be solved, in theory, by finding the fundamental matrix $(t) which is a solution to the homogeneous first order equation d+ - = Ati dt
(Reference 3, Chapter 7).
$(t) is a nonsingular matrix since it is composed of Y linearly independent solutions, .x,(t), i = 1, 2, . . ) Y, of the system. Using this and finding a particular solution of the corresponding non-homogeneous system eventually yields D(t) =
rC/(r)$-'(0)0(O) +
$(t)J'$-
V(r) = Ii/(t)j’;L- ‘(T)S(r) dr 0
0
'(7)S(7) dz
since o(O) = 0.
Numerical simulation of a compartmental
The fact that det(A) # 0 when r > 0.0 I 0 I 1 is easily
r
(1 -@aQ,+ Ml
j(
system
37
shown since
-jfi($j’
QI +QT M8
This is a direct method of finding u(t). However, finding the fundamental solution matrix for this system, although theoretically possible, would be extremely time-consuming and difficult. In addition, I+- ’ must be calculated, the integration $,I+- ‘S(7)dz and the multiplication $(t) &I+-‘(7)$ 7 )d 7 would have to be performed. This process also requires that each element of A be continuous over (0, T]. If we wish to relax continuity on all of A, requiring only piecewise continuity, we can still use this solution provided: (1) A partition 0 < t, < t2 < t, < . . . < t, < T of the interval (0,7) can be found such thatall+ i= l,..., r;j= 42 ,..., r are continuous on each subinterval (t,,tl + 1), with each of the t,‘s being a point of finite discontinuity of the aij. (2) The fundamental solution matrix can be formed for each subinterval. (3) After forming $(t) and $-l(t), i?(ct) can be found by D(t) = tj(t)$This process gives the as has been previously in terms of labor and The least formidable using matrix methods.
‘(t,)D(tJ
+ t,b(t)J’ I,- ‘(t)S(T)dT II
unique continuous solution [3] to equation (10). Unfortunately, noted, if A is not constant, this process is unrealistic to apply computation time. case-that with A and s constant-is relatively simple to solve The solution which we obtain then is:
o(t) = --A-‘?? + where R is the Y x
r
for W,, 4 + 1).
R
E-‘(f)R-‘$,
matrix composed of columns of eigenvectors of A,
1 5, (with the ordering of the eigenvalues l; corresponding in R) and
E- ‘(r) =
to the order of the eigenvectors
SHARON E. RODECAPand F. T. LINDSTROM
38
Cost and accuracy comparison of other methods to this exact method of solution when A and S are constant (using, however, standard numerical techniques for obtaining the eigenvectors and eigenvalues) will be given in Section 6. 3. NUMERICAL Many methods exist for approximating
SIMULATION the solution of matrix equations of the type
numerically. Some of the more well known are: (1) one step methods such as Euler’s method; (2) the classical Pade approximations, of which some of the best known are the backward, forward and central difference methods; (3) Taylor series and RungeKutta methods; (4) multistep predictor-corrector methods such as Milne’s method and (5) iterative methods (applied locally) such as the Gauss-Seidell or successive overrelaxation methods. Most of these methods suffer from stability problems unless we use a relatively small time step At. Since we want our simulation to be valid for long time runs (e.g. 60 days), we wish to choose a fairly large time step to keep the cost of machine time as low as possible, while maintaining good accuracy. For comparison we chose the backward difference, forward difference, central difference, Taylor series and successive over-relaxation methods. The following paragraphs describe each method. A. Backward ~$+LY~KY~ method dV ==AU+S, is approximated
as u-n+ I - (J” = A~“+i + s”+‘, At
or (1 - A At)~“”
= 0,’ + s””
At,
and u-‘I+ ] = (1 - A At)- ‘(OH + p+ I At),
(12)
if (I-AAt)- ’ exists. As will be shown later this inverse does exist and can be represented relatively easily. Also we will show that the backward difference method is unconditionally stable for matrices of the form used here. B. Forward diference method dD dt = AU + s”,
or
u-II+
’ = (I + A At@" + s” At.
Numerical simulation of a compartmental
39
system
The forward difference method is stable for 0 I Ar I
min 2rr(-‘i), ,?z,sr~
(Reference 4, p. 265).
(It will be shown later that A has negative eigenvalues, so -A has positive eigenvalues.) C. Central d@erence method
is simulated by
(13) That I_AAr
-I
( -1 2
does exist and can be represented explicitly will be proved later. The central difference method (also called the Crank-Nicholson Method) is unconditionally stable for matrices of the type used here (Reference 4, p. 265). D. Taylor series method dU -=Au+$ dr with all elements of A and S constant over [0, q, T > 0. The second derivative calculated to be d’D = A(AU + s) = A’0 + AS. dr’ In general (by induction) dkiJ ~ = A“0 + Ak-‘s drk Then we can compute 8”’ ’ by u-!I+’
=
qr,#+,)
=
D"
+
p (AjD" x
+
A'-'S")At'+
3
.i!
j-l
where p is the order of the method and the remainder R I At”+’ A”+‘M.I (J?+l)!
R"
+ A”S ’
is
SHARONE. RODECAP and F. T. LINDSTROM
40
Since 0 has a bounded steady state value and the factorial term dominates, this remainder goes to zero as p approaches infinity. Thus this series does converge for all values of t. However, for large values of At one must use many terms to large to obtain a reasonable approximation. E. Successive
over.relaxation
method (applied
Using the same approximation method,* i.e .
locally)
for du/dt
= AU + s as for the backward difference
’ = ii” + s”’ ’ At,
(I - /IA@“+
let B = I - AAt = D + E + F, where D. E and F are respectively strictly diagonal, lower and upper triangular matrices. Since this is an iterative method let r= p + ‘, K = I? + ??’ ’ ‘At and for each II we have the equation BY=
l?,
which can be changed to an iterative equation OF”’ + ‘) = -(E or
+ F)p”‘) + K (point Jacobi),
j&n + II _ - -D-
‘(E + F)?+“’ + D- ‘K.
This iterative equation can then be used to find ?= u” + ‘. To increase the speed with which y converges, we can use a relaxation factor w and obtain the following equation for each element .vi of Y:
This method is known to converge at least for w = 1, and contains the famous GaussSeidel method as a special case. (Reference 4, p. 73). It was stated without proof earlier that the matrices (I - A At) and [Z - (AA@)] are nonsingular and have inverses relatively easy to find. One can establish by direct computation that det(l - A At) and det[Z - (A At/2)] are both non-zero. Theorem 3.1 will give the actual inverse of a general class of matrices which includes both (I - A&) and [I - (A Ar/2)]. T/~o~nr 3.1. Let hII
h,.?
h 21
I722
h II 0
h 2r
h 33
B=
h 3r
.. 0
hkk
hkr -.
0
h,z
h,.,
.
.
h,,
.
. .
h,,
_
* Since it is based on the backward difference approximation, the SOR method will be limited in accuracy by the accuracy of the backward difference method. This is an unavoidable complication in adapting SOR to a system of this type.
Numerical simulation of a compartmental
41
system
The inverse of B is of the form: -
-b,,
h
/+R,H2 is -h2, RIH, h $+ RzHz T+ $+R,H,
B-l=
R3Hz
RsH,
...
R,Hl
R,H,
R,Ha
RdT
I+R,Hj
R,H,
h33
R,H,
..
W3
W,
1
kR,H,
+ Ir-
R,- ,H,- I
I
R,H,
. . .
R,H,
r24 with coeflicients determined as follows: =
9,
(e)
(a)
H,
(b)
HI!= "-+,
(c) Hjzb+,i=3
(f) R2 = ,....
r-
I,
1, (d)
H,= -;,
(j,
6 = h,,hzz
R, = bb- $RL' 11 ~,l~Zr - hlblr . -:.. /,T;;’ 1 ’
hir
(g)R.=I (h)
.,r
-’
1,
R,= -1,
- hz,b,z.
Praof The following constructive proof is based on the method of obtaining inverses given by Albert [S] which uses elementary row operations on B and the identity matrix, converting B into the identity and the identity into B -'.The following notation will be used: (i) + &)-replace the ith row by its sum with a times the jth row; a(i)-multiply the ith row by a. We start with the augmented matrix (BII): b Ir b 2r b Jr
0
b,z
br3 . . . brr
1 1
0
0 1,.
SHARON E. RODECAP and F. T. LINDSTROM
42
Performing, in order, the operations: (1)/b, ,, (2) - h2r(f)+ h11(2)/& (1) - hr2(.W7, I* and (I.) - h,,(2) with ri = h, ,h12 - h,,hz2 results in
0
b,,hzI
-h,h,
(5
6
1
0
. ..o
Then U)/h33,(4
-
~,3(3),W/‘h,,,(~)
b,.,(4),...,(~ - 1)/b,- ,r- ,,(r) - b,,- ,(r - I),
-
gives us:
R, 0 ’
1
R1
bI
I
Ci
0 1 h33
I b44
0
Rr- I
1
with Ri, i= l,..., r - 1 and y defined as in the theorem. To complete the conversion of B into the identity we simply do the following operations: (r)/y; (i) - R,(r), i = 1,. , . ,r - 1. Letting the Hi’s be as defined in the theorem, the result follows. This theorem demonstrates the existence of the needed inverse, provided bii + 0, i= I.....1-,; + 0 and h, , f b, 1. By using the explicit representation given, much computation time and some roundoff error can be eliminated from the numerical solution using either the backward difference or central difference method.
Numerical simulation of a compartmental
system
43
4. STABILITY One of the strong advantages of the backward difference method is its unconditional stability. To show this, we first need to state some definitions and develop several related theorems. Dejnition 4.1. Let A be an I x r complex matrix, and consider r distinct points or nodes, I Pi :. avery non-zero entry Ui,j of A connect Pi to P,j by means of a directed line or path Pip,. The resultant figure is called a directed graph of A. DeJinition 4.2. A directed graph is strongly connected if, for any ordered pair of nodes Pi and P, there exists a directed path . Pip,,. . . P,~.. . Pj,
connecting them. Our matrix A has the general form
r*
* *
where each of the asterisks represents
a non-zero
entry. The directed graph of A is:
This is easily seen to be strongly connected. Dejinition 4.3. An II x n (n 2 2) matrix A is reducible if there exists a permutation matrix P (i.e. it has only one nonzero entry in each row and column and each of those entries is equal to 1) such that
where A,,, is an m x m submatrix and Az,z is an (n - m) x (n - m) submatrix. If there is no such matrix A is said to be irreducible. By submatrix we mean a block consisting of the intersection of certain rows kj, j = 1,2,. . . , m with the columns of the same numbers. Theorem 4.1. An r x r matrix is irreducible if and only if its directed graph is strongly connected.
SHARON E. RODECAP and F. T. LINDSTROM
44
Proof (1) Irreducible * strongly connected. We will prove the contrapositive of this, i.e. not strongly connected * reducible. If a directed graph G(A) is not strongly connected, then we can find a maximal strongly connected subgraph F(A). Then all nodes not included in this subgraph do not have arrows both coming from and returning to the subgraph. Let (P,IF= 1 be the set of nodes in F(A). For any j,s m < j, < r write the permutation matrix Pk which will interchange each row ,jkwith a row I, 1 5 I < m3 Pl Q!F(A). For each k3 j, I m, let Pk = identity matrix. Let
This P will be the appropriate
permutation
matrix to show that A is reducible, i.e.
PAPT = (2) Strongly connected * irreducible. Again we will prove the contrapositive, ible * not strongly connected. Let A be reducible. Then 3 a permutation matrix Ps
i.e. reduc-
If A, ,i consists of rows ik, k = I,.. .,m intersected with columns .jk, k = I,..., m, then the zero elements represent the elements in rows ir, il # ib k=l , . . . , m intersected with columns ,jk, k = l,,.., m. That is, aid&= 0, k = l,.. .,m. Thus, we will have in the directed graph G(A), nodes (P,, 1 which have no arrows going to nodes [Pi,;. Since iP,,) u {Pi,) represent all the nodes, there is no possible path from a node in iPi,) to node in [P,J and therefore G(A) is not strongly connected. As a result of this theorem, we know that our matrix A is irreducible. Note also the fact that AT is irreducible follows immediately. Examining A we discover that A itself does not have diagonal dominance, but AT does, i.e.
Dtlfinition 4.4. An r x r complex matrix C is irreducibly diagonally dominant is irreducible and diagonally dominant with strict inequality in at least one i. .. Since in row 2 of AT we have a,z
T
=
@I + (21+~.QHI + r) > M2
QWQHl ) Q, + Q2 MZ
since f > 0, AT is irreducibly diagonally dominant. Theorem 4.2. Let C be an r x r complex matrix and let
ML
’
if C
Numerical
simulation
of a compartmental
system
45
then all eigenvalues of C will lie in the union of disks 1: - cii 1 I Ai, I I i I r. Since this is a standard theorem we will omit the proof. The following is a result from Varga [4] which sharpens the Gerschgorin theorem. Theorrm 4.3. Let C be an irreducible r x Y complex matrix. If A, an eigenvalue of C, is a boundary point of the union of disks 1z - Cii 1 I Ai, then all r circles 1z - Cii I = /li pass through i. Proqf Let A be an eigenvalue of C and _Yits corresponding eigenvector, normalized so that 1.~~1= 1 2 1.~~1,for all 1 I i 5 r where .yk is the largest component in absolute value. Then
If 1.~~1 = 1 we have
If 2 is assumed to be a boundary
point of the union of all the disks 1,”-
CjiI I
Ai,
then we must have
and. thus
This requires that 1.~~1= 1 for all ckj # 0. Since C is irreducible, for at least one p # k, ckp # 0; and hence, J.Y,,/= 1. Repeating the above argument with k = p we obtain
1~- cppl = ~ ICpjJ.
j- I
/.\.,jl
=
np.
i+r,
again using the hypothesis that 1 is a boundary point of the union of all disks IZ - Cii( I ,li. Therefore, 1.~~1 = 1V C,j # 0. By Theorem 4.1 since C is strongly connected we can find a sequence of non-zero entries of C such as C,,, cpk,,. . . , ck,_,, for any ,j, 1 < ,j I I’. This allows us to continue our argument and conclude that 11 - czjjI = Aj for all I I ,j I I’,
demonstrating that all r circles 1: - cjil = “j pass through A. And, finally. another theorem from Varga [4] gives more about the nature of the 2-S.
Theorem 4.4. Let C be an r x r irreducibly diagonally dominant complex matrix. If the diagonal entries of C are all positive real numbers, then for all the eigenvalues of C, re(lJ > 0, 1 I i I r.
46
Pro@‘: By the Gerschgorin
SHARON E. RODECAPand F. T. LINDSTROM
theorem all eigenvalues of C lie in the union of disks,
Since C is irreducibly diagonally dominant, Ic!&l2 A,Vk and Icu I > A, for some 1. Let 2 be an eigenvalue of C. Then 1 is in at least one disk D, = [.zllz - ckkl I A,). Suppose re(2) = 0, and 1 = hi.h # 0. (hi - ckkl = , czi y p > ,I, and & and 1,# L; D,, a contradiction. If h = 0, then IA - ckk(= (cJ 2 /1, so since /1 E D,, /cJ = A, and i = 0 is a boundary point of D,. By Theorem 4.3, 0 is then a boundary point of all disks D,, k = I,..., r. However. 0 cannot be a boundary point of D, since (0 - cIll = lqlJ > A,. Therefore, 0 cannot be an eigenvalue of C (also showing C to be nonsingular). If re(l) < 0, i = a + b, then 13,- ckkl = &
- cJ2 + b2 > /i, V k,
since a < 0, ckk > 0; and again, A$ u D,, a contradiction.
Thus re(2) > 0 for all eigenvalues of C. Neither A nor AT satisfy the criterion in Theorem 4.4. However, -(AT) = (-AT) does. Thus, as the eigenvalues of (-AT) and -A are the same (Ref. 6, p. 93), we know that -A also has eigenvalues with positive real parts by Theorem 4.4. Dclfinition 4.5. A matrix approximation of dU -~ = Ali + s, dt
(10)
in the form of u II + 1 = T(b”’ + s” + ‘At) is unconditionally stable if p(T) < 1 (where p(T) is the spectral radius of T) for all t > 0 (i.e. it is stable for any At). Theorem 4.5 [4.]. Let C be an r x r matrix whose eigenvalues all have positive real parts. Then the matrix approximation T= (I + CAt)- ’ for the backward difference method is unconditionally stable. ProuJ: If the eigenvalues of C are I&:; , then the eigenvalues of (I + CAt) are I I + iiAt). Re(1 + li At) > 1 Vi. Similarly the eigenvalues of (I + CAt)- ’ are [ 1 + ,$At) - ‘. Since re (1 + AiAt) > 1, p(l + CAt) > 1 and &(I + CAt)- ‘1 < 1 which insures that the backward difference method is unconditionally stable when the real parts of the eigenvalues of C are positive. If we let C = -A in Theorem 4.5, we see that our matrix approximation with (I AAt)- ’ is unconditionally stable. 5. EXAMPLE
As an example of the use of the system we have developed here, the following case was run. The mammal used is a typical 300g mature Wistar@ rat, previously fed on
Numerical
simulation
of a compartmental
system
47
a dieldrin-free Purina Chow@ diet. We assume that all Q’s and M’s are fixed in time and that the dieldrin destruction rate coefficient (r) is constant. In addition, suppose that the rat is eating about 20 g of food daily [7] over a 12 hr period and the dieldrin concentration in the diet is 25 pg dieldrin/g food; that is, 9 = 41.67 pg/hr over the I2-hr feeding period and 9 = 0 for the remaining 12 hr in which the rat essentially eats nothing. Thus Y is a periodic function of time which can be extended for as many days as is desired. The tissue lipid fractions V;) used in calculation of the compartmental lipid masses (total mass.i, = M,) are derived from Long [9], Williams and Landsford [9] and Spector [lo]. This calculation is demonstrated in Table 1 along with the flow rates of blood lipid in each of the compartments (Qs). After cu. 4 wk on this feeding and dosage schedule, the adipose tissue dieldrin concentrations build up in an oscillatory fashion to almost a “steady state”. Figure 3 shows a plot of the lipid phase dieldrin concentration in blood and adipose tissues generated by the values in Table 1, when the rat is fed on the schedule mentioned above for 40 days to “steady-state” conditions and then put on a dieldrin-free diet for 20 more days. In order to generate these curves, an estimate of r (constant in this example) must be obtained, as explained below. At “steady state” conditions, after feeding the rats for several weeks, Walker et al. [ 1l] and Barron and Walton [ 121 found a depot fat (whole tissue) dieldrin concentration of about 50 pg dieldrin/g adipose tissue at a 25 ppm dietary concentration of dieldrin. If we assume a lipid fraction fF of 0.637 (g lipid/gadipose tissue), and we know that the depot fat lipid phase dieldrin concentration CIF (pg dieldrin/g lipid) is related to the depot fat tissue dieldrin concentration CF (pg dieldrin/g adipose tissue) by
Table
1. Compartment
mass, blood
flow and other important 300 g mature male rat
parameters
for a hypothetical
Blood Fraction
Ltpld bias* x
Symbol
compartment
7.
Blood pool B (blood chambers of heart. lungs. arteries, and veins1
f
(gms)
ms tg gms
23.4
Lipid M
x
lip*d
tissue)
kms
mow
Mass
4x
x
lipid)
0.004
0.09
Rate
(
gms
hrblood)
Blood
Lipid
Flow
,gms
Rate
Qx blood hr
4,751.b
18.21
Ktdney
K
2.7
0 02
0.05
698.0
2.79
G~St~0*llt~S!l"~l (small bow=!
G.I.
15,s
0.07
1. 11
811.1
3.24
3.96
LiVer
I.
11 9
0 06
0 72
989.3
Bra,"
BR
1.9
0.09
0. 17
680. I
2.72
DepJt Fat
F
54.0
0. 64
34.56
710.1
2.84
Muscle (includes associated fat)
M
I20.0
0.07
8.40
533.3
2.13
6 60
133.3
0.53
Other tis8"e OT 60.0 0. 11 ,sk,n, etc., Rex = 0.0 (Heath and Vandakar. 19641 = 100.0 (estimated from data of Cole et al. 1197011 B
e e
= =
0.995. 0.050,
’ Lt < 24hours I2 hours1 repeatable
12
‘_t
daily
(periodtc function)
0. For these particular flow and mass parameter values a short summary of the accuracy of the numerical simulation schemes is presented in Table 2. Each method is compared to the “exact” method which was stated in Section 2 (using the standard inverse and eigenvalue-eigenvector routines of the OSU Computer Center). lPUC _
Al
Backward Central
Forward
20 hrs
,218
.0539
difference
,595
difference
SOR Taylor
order
= 10.0hr
10 klre
IO hrs
0318
.0539
““ataMe
far
20 hrs
0170
3457 .23X3
.,218
1 Ohr
dt= 30 hrs
-0318
’
lPC’ I
I_
Biol. Med.
Pergamon
Press 1976. Vol. 6, pp. 33-51.
Printed in Great
Britain.
NUMERICAL SIMULATION OF A COMPARTMENTAL SYSTEM FOR THE DISTRIBUTION OF LIPID SOLUBLE CHEMICALS IN MAMMALIAN TISSUES SHARON E. RODECAP Departments
and F. T. LINDSTROM
of Mathematics and Statistics, Oregon State University, Corvallis, Oregon 97331, U.S.A. (Received 23 October 1974)
Abstract-An analysis of the stability and accuracy of a particular numerical scheme used in the mathematical modeling of the flow limit pharmacodynamics of highly lipid soluble drugs in mammalian tissue systems is presented. The necessary definitions and theorems together with their proofs are included. A typical example of the flow limit pharmacodynamics of the highly lipid soluble drug dieldrin (HEOD) is also presented. Numerical dieldrin
Mathematical model
Flow limit pharmacodynamics
Lipid soluble drug
1. INTRODUCTION In recent years, with the development of high speed digital computers, the use of mathematical compartmental models to describe and predict biological phenomena has become more prevalent. One of the major uses has been to describe blood flow and drug distribution in mammalian tissue systems. A compartmental model is set up by partitioning the body into a finite number of compartments [ 11, each of which consists of tissues which are relatively homogeneous with respect to the drug being studied. To each compartment is assigned a variable denoting the amount of drug in that compartment. A system of differential equations is then developed giving the time rate of change of the amount of drug in each compartment in terms of the current amounts in each of the compartments, input from outside sources and excretion and destruction of the chemical. Correct simulation of these systems involves knowing a priori the compartment masses (volumes) and flow rates of the blood (in which the drug is assumed to be dissolved) between compartments. Using this basic method Lindstrom et al. [2] have formulated a preliminary model to describe the pharmacokinetics and the pharmacodynamics (drug distribution mechanics) of highly lipid (fat) soluble drugs in mammalian tissue systems. The following assumptions are made : (I ) The drug is transported via the blood stream completely dissolved in blood lipids (fats). (2) Drug transfer from the blood capillaries of the lipid fraction of each compartment is rapid (flow limit conditions operating). (3) The drug is converted (metabolised) only in the liver and is itself not easily excreted by the kidneys. 33
34
SHARON E. RODECAP and F. T. LINDSTROM
(4) Tissues can be partitioned into relatively homogeneous compartments within which there is uniform, instantaneous mixing (no intracompartment concentration gradients, difTusion currents, exist). (5) There exists a complete system of mass balance equations describing the drug distribution. (6) All blood lipid flow rates (QJ and compartmental lipid masses (M,) are continuous functions of time (0 I t I x). Note that these assumptions and the following mathematics are general enough to be useful for simulation of any highly lipid soluble drug in the body. Figure 1 is a diagram of the system with rectangles representing compartments, I/;
I
’
(Fecal excretion)
Fig. 1.
the amount of drug in compartment i (g), Mi the lipid mass (g), Qi the flow rate of blood lipid through the ith compartment (g/hr), 9 the input (source) function (g/hr). and r the liver drug conversion rate coefficient (g/hr). The remaining letters represent various biological constants. The lumen, although shown as a rectangle, is not a regular compartment; but it is simply a part of the liver compartment into which the drug enters with some slight excretion of the drug allowed. The term “mass balance equation” in assumption (5) simply means h:
the time rate of change of drug mass dnJdt i
iI
= [drug mass flow delivered via arteries from the general blood pool] + [drug sources] - [drug conversion, (sinks)] E) ( - [drug mass removed from compartment via venous blood]
Numerical simulation of a compartmental
system
Using this basic concept of mass balance on each compartment, equations is derived:
&L dt
0 .y + S?QH,f$ z
du,
z=~,z+Qs
-
8 du, dt = Q32
!+
Q,
2 -
Q,
35
the following system of
nil;-, (G.I.),
x
(aQHI + Q, + Qz + r)$
- (Q3 + R,,)$
(1)
I
Z
(Liver),
(Kidney),
(2) (3)
3
$f
=
Q4s
x
- Q3$
Q5$
- Q,$
du, dt - Qh 2
- Qh$
2
=
du, dt = Q7 $
x
(Brain),
(4)
(Depot fat),
(5)
(Muscle),
(6)
(Other tissues),
(7)
4
- Q, $
du, dr - (Q, + Q,,$
2
5 h 7
+ i ai .j=3
- Q8$
(Blood pool), 8
I
(8)
with the auxiliary equation: 7
Q8 = z
Qj (Conservation
of
flow).
j-l
(9)
Equations (l-8) can be written - = AD + s(t). dt
(10)
where OF-
I:_ 0
0
S=
.
9
0
and the matrix A appears in Fig. 2. Note that all the mass and flow parameters Mi and (2; respectively, are assumed to be real positive piecewise continuous functions of time t over CO,7J, T> 0.The special terms, 13, R,,,, 9, and r are positive real parameters with 0 I 8 I 1. The term 9. the source function, is assumed to be at least bounded and measurable [0, 13. For the most part R,,,is simply assumed equal to zero.
SHARON
E.
RODECAP
and F. T.
LINDS~OM
-Q1
-5 Ql
T -(Q3+Rexl 0
0
0
M3 0
0
0
0
0
0
0
0
0
Ql+Qz --
-Q4 M,
0
0
0
0
0
0
0
“b M6
0
0
0
0
-Q5 M5
0
0
0
0
M2
Fig. 2. Mammalian transport matrix for a male mammal.
2.
FORMAL
SOLUTION
Consider the first order system dU - = AD + S(r), 0 I t < cx, dt D(O) = 0, det(A) # 0,
with a;,i = aij(t), i= 1,2,. . . ,r; j = 1,2,. . . , Y continuous [O 5 t < Y&T > 0) and with all elements of S piecewise continuous on the same interval. This can be solved, in theory, by finding the fundamental matrix $(t) which is a solution to the homogeneous first order equation d+ - = Ati dt
(Reference 3, Chapter 7).
$(t) is a nonsingular matrix since it is composed of Y linearly independent solutions, .x,(t), i = 1, 2, . . ) Y, of the system. Using this and finding a particular solution of the corresponding non-homogeneous system eventually yields D(t) =
rC/(r)$-'(0)0(O) +
$(t)J'$-
V(r) = Ii/(t)j’;L- ‘(T)S(r) dr 0
0
'(7)S(7) dz
since o(O) = 0.
Numerical simulation of a compartmental
The fact that det(A) # 0 when r > 0.0 I 0 I 1 is easily
r
(1 -@aQ,+ Ml
j(
system
37
shown since
-jfi($j’
QI +QT M8
This is a direct method of finding u(t). However, finding the fundamental solution matrix for this system, although theoretically possible, would be extremely time-consuming and difficult. In addition, I+- ’ must be calculated, the integration $,I+- ‘S(7)dz and the multiplication $(t) &I+-‘(7)$ 7 )d 7 would have to be performed. This process also requires that each element of A be continuous over (0, T]. If we wish to relax continuity on all of A, requiring only piecewise continuity, we can still use this solution provided: (1) A partition 0 < t, < t2 < t, < . . . < t, < T of the interval (0,7) can be found such thatall+ i= l,..., r;j= 42 ,..., r are continuous on each subinterval (t,,tl + 1), with each of the t,‘s being a point of finite discontinuity of the aij. (2) The fundamental solution matrix can be formed for each subinterval. (3) After forming $(t) and $-l(t), i?(ct) can be found by D(t) = tj(t)$This process gives the as has been previously in terms of labor and The least formidable using matrix methods.
‘(t,)D(tJ
+ t,b(t)J’ I,- ‘(t)S(T)dT II
unique continuous solution [3] to equation (10). Unfortunately, noted, if A is not constant, this process is unrealistic to apply computation time. case-that with A and s constant-is relatively simple to solve The solution which we obtain then is:
o(t) = --A-‘?? + where R is the Y x
r
for W,, 4 + 1).
R
E-‘(f)R-‘$,
matrix composed of columns of eigenvectors of A,
1 5, (with the ordering of the eigenvalues l; corresponding in R) and
E- ‘(r) =
to the order of the eigenvectors
SHARON E. RODECAPand F. T. LINDSTROM
38
Cost and accuracy comparison of other methods to this exact method of solution when A and S are constant (using, however, standard numerical techniques for obtaining the eigenvectors and eigenvalues) will be given in Section 6. 3. NUMERICAL Many methods exist for approximating
SIMULATION the solution of matrix equations of the type
numerically. Some of the more well known are: (1) one step methods such as Euler’s method; (2) the classical Pade approximations, of which some of the best known are the backward, forward and central difference methods; (3) Taylor series and RungeKutta methods; (4) multistep predictor-corrector methods such as Milne’s method and (5) iterative methods (applied locally) such as the Gauss-Seidell or successive overrelaxation methods. Most of these methods suffer from stability problems unless we use a relatively small time step At. Since we want our simulation to be valid for long time runs (e.g. 60 days), we wish to choose a fairly large time step to keep the cost of machine time as low as possible, while maintaining good accuracy. For comparison we chose the backward difference, forward difference, central difference, Taylor series and successive over-relaxation methods. The following paragraphs describe each method. A. Backward ~$+LY~KY~ method dV ==AU+S, is approximated
as u-n+ I - (J” = A~“+i + s”+‘, At
or (1 - A At)~“”
= 0,’ + s””
At,
and u-‘I+ ] = (1 - A At)- ‘(OH + p+ I At),
(12)
if (I-AAt)- ’ exists. As will be shown later this inverse does exist and can be represented relatively easily. Also we will show that the backward difference method is unconditionally stable for matrices of the form used here. B. Forward diference method dD dt = AU + s”,
or
u-II+
’ = (I + A At@" + s” At.
Numerical simulation of a compartmental
39
system
The forward difference method is stable for 0 I Ar I
min 2rr(-‘i), ,?z,sr~
(Reference 4, p. 265).
(It will be shown later that A has negative eigenvalues, so -A has positive eigenvalues.) C. Central d@erence method
is simulated by
(13) That I_AAr
-I
( -1 2
does exist and can be represented explicitly will be proved later. The central difference method (also called the Crank-Nicholson Method) is unconditionally stable for matrices of the type used here (Reference 4, p. 265). D. Taylor series method dU -=Au+$ dr with all elements of A and S constant over [0, q, T > 0. The second derivative calculated to be d’D = A(AU + s) = A’0 + AS. dr’ In general (by induction) dkiJ ~ = A“0 + Ak-‘s drk Then we can compute 8”’ ’ by u-!I+’
=
qr,#+,)
=
D"
+
p (AjD" x
+
A'-'S")At'+
3
.i!
j-l
where p is the order of the method and the remainder R I At”+’ A”+‘M.I (J?+l)!
R"
+ A”S ’
is
SHARONE. RODECAP and F. T. LINDSTROM
40
Since 0 has a bounded steady state value and the factorial term dominates, this remainder goes to zero as p approaches infinity. Thus this series does converge for all values of t. However, for large values of At one must use many terms to large to obtain a reasonable approximation. E. Successive
over.relaxation
method (applied
Using the same approximation method,* i.e .
locally)
for du/dt
= AU + s as for the backward difference
’ = ii” + s”’ ’ At,
(I - /IA@“+
let B = I - AAt = D + E + F, where D. E and F are respectively strictly diagonal, lower and upper triangular matrices. Since this is an iterative method let r= p + ‘, K = I? + ??’ ’ ‘At and for each II we have the equation BY=
l?,
which can be changed to an iterative equation OF”’ + ‘) = -(E or
+ F)p”‘) + K (point Jacobi),
j&n + II _ - -D-
‘(E + F)?+“’ + D- ‘K.
This iterative equation can then be used to find ?= u” + ‘. To increase the speed with which y converges, we can use a relaxation factor w and obtain the following equation for each element .vi of Y:
This method is known to converge at least for w = 1, and contains the famous GaussSeidel method as a special case. (Reference 4, p. 73). It was stated without proof earlier that the matrices (I - A At) and [Z - (AA@)] are nonsingular and have inverses relatively easy to find. One can establish by direct computation that det(l - A At) and det[Z - (A At/2)] are both non-zero. Theorem 3.1 will give the actual inverse of a general class of matrices which includes both (I - A&) and [I - (A Ar/2)]. T/~o~nr 3.1. Let hII
h,.?
h 21
I722
h II 0
h 2r
h 33
B=
h 3r
.. 0
hkk
hkr -.
0
h,z
h,.,
.
.
h,,
.
. .
h,,
_
* Since it is based on the backward difference approximation, the SOR method will be limited in accuracy by the accuracy of the backward difference method. This is an unavoidable complication in adapting SOR to a system of this type.
Numerical simulation of a compartmental
41
system
The inverse of B is of the form: -
-b,,
h
/+R,H2 is -h2, RIH, h $+ RzHz T+ $+R,H,
B-l=
R3Hz
RsH,
...
R,Hl
R,H,
R,Ha
RdT
I+R,Hj
R,H,
h33
R,H,
..
W3
W,
1
kR,H,
+ Ir-
R,- ,H,- I
I
R,H,
. . .
R,H,
r24 with coeflicients determined as follows: =
9,
(e)
(a)
H,
(b)
HI!= "-+,
(c) Hjzb+,i=3
(f) R2 = ,....
r-
I,
1, (d)
H,= -;,
(j,
6 = h,,hzz
R, = bb- $RL' 11 ~,l~Zr - hlblr . -:.. /,T;;’ 1 ’
hir
(g)R.=I (h)
.,r
-’
1,
R,= -1,
- hz,b,z.
Praof The following constructive proof is based on the method of obtaining inverses given by Albert [S] which uses elementary row operations on B and the identity matrix, converting B into the identity and the identity into B -'.The following notation will be used: (i) + &)-replace the ith row by its sum with a times the jth row; a(i)-multiply the ith row by a. We start with the augmented matrix (BII): b Ir b 2r b Jr
0
b,z
br3 . . . brr
1 1
0
0 1,.
SHARON E. RODECAP and F. T. LINDSTROM
42
Performing, in order, the operations: (1)/b, ,, (2) - h2r(f)+ h11(2)/& (1) - hr2(.W7, I* and (I.) - h,,(2) with ri = h, ,h12 - h,,hz2 results in
0
b,,hzI
-h,h,
(5
6
1
0
. ..o
Then U)/h33,(4
-
~,3(3),W/‘h,,,(~)
b,.,(4),...,(~ - 1)/b,- ,r- ,,(r) - b,,- ,(r - I),
-
gives us:
R, 0 ’
1
R1
bI
I
Ci
0 1 h33
I b44
0
Rr- I
1
with Ri, i= l,..., r - 1 and y defined as in the theorem. To complete the conversion of B into the identity we simply do the following operations: (r)/y; (i) - R,(r), i = 1,. , . ,r - 1. Letting the Hi’s be as defined in the theorem, the result follows. This theorem demonstrates the existence of the needed inverse, provided bii + 0, i= I.....1-,; + 0 and h, , f b, 1. By using the explicit representation given, much computation time and some roundoff error can be eliminated from the numerical solution using either the backward difference or central difference method.
Numerical simulation of a compartmental
system
43
4. STABILITY One of the strong advantages of the backward difference method is its unconditional stability. To show this, we first need to state some definitions and develop several related theorems. Dejnition 4.1. Let A be an I x r complex matrix, and consider r distinct points or nodes, I Pi :. avery non-zero entry Ui,j of A connect Pi to P,j by means of a directed line or path Pip,. The resultant figure is called a directed graph of A. DeJinition 4.2. A directed graph is strongly connected if, for any ordered pair of nodes Pi and P, there exists a directed path . Pip,,. . . P,~.. . Pj,
connecting them. Our matrix A has the general form
r*
* *
where each of the asterisks represents
a non-zero
entry. The directed graph of A is:
This is easily seen to be strongly connected. Dejinition 4.3. An II x n (n 2 2) matrix A is reducible if there exists a permutation matrix P (i.e. it has only one nonzero entry in each row and column and each of those entries is equal to 1) such that
where A,,, is an m x m submatrix and Az,z is an (n - m) x (n - m) submatrix. If there is no such matrix A is said to be irreducible. By submatrix we mean a block consisting of the intersection of certain rows kj, j = 1,2,. . . , m with the columns of the same numbers. Theorem 4.1. An r x r matrix is irreducible if and only if its directed graph is strongly connected.
SHARON E. RODECAP and F. T. LINDSTROM
44
Proof (1) Irreducible * strongly connected. We will prove the contrapositive of this, i.e. not strongly connected * reducible. If a directed graph G(A) is not strongly connected, then we can find a maximal strongly connected subgraph F(A). Then all nodes not included in this subgraph do not have arrows both coming from and returning to the subgraph. Let (P,IF= 1 be the set of nodes in F(A). For any j,s m < j, < r write the permutation matrix Pk which will interchange each row ,jkwith a row I, 1 5 I < m3 Pl Q!F(A). For each k3 j, I m, let Pk = identity matrix. Let
This P will be the appropriate
permutation
matrix to show that A is reducible, i.e.
PAPT = (2) Strongly connected * irreducible. Again we will prove the contrapositive, ible * not strongly connected. Let A be reducible. Then 3 a permutation matrix Ps
i.e. reduc-
If A, ,i consists of rows ik, k = I,.. .,m intersected with columns .jk, k = I,..., m, then the zero elements represent the elements in rows ir, il # ib k=l , . . . , m intersected with columns ,jk, k = l,,.., m. That is, aid&= 0, k = l,.. .,m. Thus, we will have in the directed graph G(A), nodes (P,, 1 which have no arrows going to nodes [Pi,;. Since iP,,) u {Pi,) represent all the nodes, there is no possible path from a node in iPi,) to node in [P,J and therefore G(A) is not strongly connected. As a result of this theorem, we know that our matrix A is irreducible. Note also the fact that AT is irreducible follows immediately. Examining A we discover that A itself does not have diagonal dominance, but AT does, i.e.
Dtlfinition 4.4. An r x r complex matrix C is irreducibly diagonally dominant is irreducible and diagonally dominant with strict inequality in at least one i. .. Since in row 2 of AT we have a,z
T
=
@I + (21+~.QHI + r) > M2
QWQHl ) Q, + Q2 MZ
since f > 0, AT is irreducibly diagonally dominant. Theorem 4.2. Let C be an r x r complex matrix and let
ML
’
if C
Numerical
simulation
of a compartmental
system
45
then all eigenvalues of C will lie in the union of disks 1: - cii 1 I Ai, I I i I r. Since this is a standard theorem we will omit the proof. The following is a result from Varga [4] which sharpens the Gerschgorin theorem. Theorrm 4.3. Let C be an irreducible r x Y complex matrix. If A, an eigenvalue of C, is a boundary point of the union of disks 1z - Cii 1 I Ai, then all r circles 1z - Cii I = /li pass through i. Proqf Let A be an eigenvalue of C and _Yits corresponding eigenvector, normalized so that 1.~~1= 1 2 1.~~1,for all 1 I i 5 r where .yk is the largest component in absolute value. Then
If 1.~~1 = 1 we have
If 2 is assumed to be a boundary
point of the union of all the disks 1,”-
CjiI I
Ai,
then we must have
and. thus
This requires that 1.~~1= 1 for all ckj # 0. Since C is irreducible, for at least one p # k, ckp # 0; and hence, J.Y,,/= 1. Repeating the above argument with k = p we obtain
1~- cppl = ~ ICpjJ.
j- I
/.\.,jl
=
np.
i+r,
again using the hypothesis that 1 is a boundary point of the union of all disks IZ - Cii( I ,li. Therefore, 1.~~1 = 1V C,j # 0. By Theorem 4.1 since C is strongly connected we can find a sequence of non-zero entries of C such as C,,, cpk,,. . . , ck,_,, for any ,j, 1 < ,j I I’. This allows us to continue our argument and conclude that 11 - czjjI = Aj for all I I ,j I I’,
demonstrating that all r circles 1: - cjil = “j pass through A. And, finally. another theorem from Varga [4] gives more about the nature of the 2-S.
Theorem 4.4. Let C be an r x r irreducibly diagonally dominant complex matrix. If the diagonal entries of C are all positive real numbers, then for all the eigenvalues of C, re(lJ > 0, 1 I i I r.
46
Pro@‘: By the Gerschgorin
SHARON E. RODECAPand F. T. LINDSTROM
theorem all eigenvalues of C lie in the union of disks,
Since C is irreducibly diagonally dominant, Ic!&l2 A,Vk and Icu I > A, for some 1. Let 2 be an eigenvalue of C. Then 1 is in at least one disk D, = [.zllz - ckkl I A,). Suppose re(2) = 0, and 1 = hi.h # 0. (hi - ckkl = , czi y p > ,I, and & and 1,# L; D,, a contradiction. If h = 0, then IA - ckk(= (cJ 2 /1, so since /1 E D,, /cJ = A, and i = 0 is a boundary point of D,. By Theorem 4.3, 0 is then a boundary point of all disks D,, k = I,..., r. However. 0 cannot be a boundary point of D, since (0 - cIll = lqlJ > A,. Therefore, 0 cannot be an eigenvalue of C (also showing C to be nonsingular). If re(l) < 0, i = a + b, then 13,- ckkl = &
- cJ2 + b2 > /i, V k,
since a < 0, ckk > 0; and again, A$ u D,, a contradiction.
Thus re(2) > 0 for all eigenvalues of C. Neither A nor AT satisfy the criterion in Theorem 4.4. However, -(AT) = (-AT) does. Thus, as the eigenvalues of (-AT) and -A are the same (Ref. 6, p. 93), we know that -A also has eigenvalues with positive real parts by Theorem 4.4. Dclfinition 4.5. A matrix approximation of dU -~ = Ali + s, dt
(10)
in the form of u II + 1 = T(b”’ + s” + ‘At) is unconditionally stable if p(T) < 1 (where p(T) is the spectral radius of T) for all t > 0 (i.e. it is stable for any At). Theorem 4.5 [4.]. Let C be an r x r matrix whose eigenvalues all have positive real parts. Then the matrix approximation T= (I + CAt)- ’ for the backward difference method is unconditionally stable. ProuJ: If the eigenvalues of C are I&:; , then the eigenvalues of (I + CAt) are I I + iiAt). Re(1 + li At) > 1 Vi. Similarly the eigenvalues of (I + CAt)- ’ are [ 1 + ,$At) - ‘. Since re (1 + AiAt) > 1, p(l + CAt) > 1 and &(I + CAt)- ‘1 < 1 which insures that the backward difference method is unconditionally stable when the real parts of the eigenvalues of C are positive. If we let C = -A in Theorem 4.5, we see that our matrix approximation with (I AAt)- ’ is unconditionally stable. 5. EXAMPLE
As an example of the use of the system we have developed here, the following case was run. The mammal used is a typical 300g mature Wistar@ rat, previously fed on
Numerical
simulation
of a compartmental
system
47
a dieldrin-free Purina Chow@ diet. We assume that all Q’s and M’s are fixed in time and that the dieldrin destruction rate coefficient (r) is constant. In addition, suppose that the rat is eating about 20 g of food daily [7] over a 12 hr period and the dieldrin concentration in the diet is 25 pg dieldrin/g food; that is, 9 = 41.67 pg/hr over the I2-hr feeding period and 9 = 0 for the remaining 12 hr in which the rat essentially eats nothing. Thus Y is a periodic function of time which can be extended for as many days as is desired. The tissue lipid fractions V;) used in calculation of the compartmental lipid masses (total mass.i, = M,) are derived from Long [9], Williams and Landsford [9] and Spector [lo]. This calculation is demonstrated in Table 1 along with the flow rates of blood lipid in each of the compartments (Qs). After cu. 4 wk on this feeding and dosage schedule, the adipose tissue dieldrin concentrations build up in an oscillatory fashion to almost a “steady state”. Figure 3 shows a plot of the lipid phase dieldrin concentration in blood and adipose tissues generated by the values in Table 1, when the rat is fed on the schedule mentioned above for 40 days to “steady-state” conditions and then put on a dieldrin-free diet for 20 more days. In order to generate these curves, an estimate of r (constant in this example) must be obtained, as explained below. At “steady state” conditions, after feeding the rats for several weeks, Walker et al. [ 1l] and Barron and Walton [ 121 found a depot fat (whole tissue) dieldrin concentration of about 50 pg dieldrin/g adipose tissue at a 25 ppm dietary concentration of dieldrin. If we assume a lipid fraction fF of 0.637 (g lipid/gadipose tissue), and we know that the depot fat lipid phase dieldrin concentration CIF (pg dieldrin/g lipid) is related to the depot fat tissue dieldrin concentration CF (pg dieldrin/g adipose tissue) by
Table
1. Compartment
mass, blood
flow and other important 300 g mature male rat
parameters
for a hypothetical
Blood Fraction
Ltpld bias* x
Symbol
compartment
7.
Blood pool B (blood chambers of heart. lungs. arteries, and veins1
f
(gms)
ms tg gms
23.4
Lipid M
x
lip*d
tissue)
kms
mow
Mass
4x
x
lipid)
0.004
0.09
Rate
(
gms
hrblood)
Blood
Lipid
Flow
,gms
Rate
Qx blood hr
4,751.b
18.21
Ktdney
K
2.7
0 02
0.05
698.0
2.79
G~St~0*llt~S!l"~l (small bow=!
G.I.
15,s
0.07
1. 11
811.1
3.24
3.96
LiVer
I.
11 9
0 06
0 72
989.3
Bra,"
BR
1.9
0.09
0. 17
680. I
2.72
DepJt Fat
F
54.0
0. 64
34.56
710.1
2.84
Muscle (includes associated fat)
M
I20.0
0.07
8.40
533.3
2.13
6 60
133.3
0.53
Other tis8"e OT 60.0 0. 11 ,sk,n, etc., Rex = 0.0 (Heath and Vandakar. 19641 = 100.0 (estimated from data of Cole et al. 1197011 B
e e
= =
0.995. 0.050,
’ Lt < 24hours I2 hours1 repeatable
12
‘_t
daily
(periodtc function)
0. For these particular flow and mass parameter values a short summary of the accuracy of the numerical simulation schemes is presented in Table 2. Each method is compared to the “exact” method which was stated in Section 2 (using the standard inverse and eigenvalue-eigenvector routines of the OSU Computer Center). lPUC _
Al
Backward Central
Forward
20 hrs
,218
.0539
difference
,595
difference
SOR Taylor
order
= 10.0hr
10 klre
IO hrs
0318
.0539
““ataMe
far
20 hrs
0170
3457 .23X3
.,218
1 Ohr
dt= 30 hrs
-0318
’
lPC’ I
I_
