Computer simulation of ferrokinetic models

COMPUTER PROGRAMS IN BIOMEDICINE 1 (1970) 90-104. NORTH-HOLLAND PUBLISHING COMPANY

COMPUTER

SIMULATION

OF FERROKINETIC

MODELS*

Torgny GROTH and Werner SCHNEIDER

UppsalaData Center, University of Uppsala, Sweden Erik SANDEWALL

Department of Computer Sciences, University of Uppsala, Sweden and Jean-Claude VUILLE

Department of Pediatrics, University Hospital University of Uppsala, Sweden

Camputer methods for analysis of iron kinetics on the basis of specified models are presented. One method is basec on direct simulation of the processes in the model, whereas the other method uses a standard fourth-order Runge - Kutta integration of the system of differential equations describing the model. Simulation

1. Introduction and computational methods Observation of the behaviour in plasma and red cells of intravenously injected radioiron form a basis for the study o f iron-kinetics in man. The theoretical interpretation of the experimental data can be derived from different ferrokinetic models consisting o f various iron pools and flows connected to the plasma compartment in various manners. In the original paper by Garby et al. [ 1] a relatively simple model was presented which could be formulated mathematically in terms of analytically solvable differential equations. This model was in many respects similar to that proposed in 1961 by Pollycove and Mortimer [2]. Further investigations [3,4], however, especially on the bone marrow part of the model, demonstrated the necessity o f analysing models o f different basic structure. The introduction o f more sophisticated mechanisms, simulating maturation o f and iron uptake by red cell precursors in the bone marrow, as well as premature death of newly formed red cells, necessitated * This work was supported in part by the Swedish Natural Science Research Council.

Ferrokinetic Models

the use o f new computer techniques, since it was no longer possible to solve the differential equations analytically. The purpose of this paper is to illustrate two suitable computer methods for the simulation of this kind o f process, to be used in the cyclic search for adequate models, including experiments, data reduction, formulation o f models, parameter estimation and simulation o f new experiments. In the first method (A) the simulation is performed in small time intervals At and the change Ax in each compartment x is calculated from the system of differential equations describing the system by approximating the derivatives by quotients o f differences Ax/At. The other method (B) is based on the solution o f the system of differential equations by a standard fourth order Runge - Kutta algorithm. For a detailde review of different methods for the simulation of continuous systems see ref. [5]. Method A has two interdependent drawbacks for routine use. In order to obtain high accuracy, the time intervals used in the simulation procedure must be very small, a requirement which increases the computer time needed for the production o f one set o f

91

COMPUTER SIMULATIONOF FERROKINETIC MODELS theoretical curves. Using method B larger time intervals can be chosen for the same level of accuracy, thereby reducing the total number of steps at which calculations have to be performed. Both methods are able to yield valid theoretical curves for a specified model, but it may be recommended to use the simple procedure of method A when several models are to be studied, e.g. with respect to their compatibility with experimental data. Once a specified model is accepted for further work, which often implies the production of a large number of theoretical curves (parameter estimation by curve fitting), a more powerful procedure like method B should be applied. From the point of view of the ease and the flexibility of the programming technique, however, method A is clearly superior to method B. Available general programs for the analysis of multi-compartment systems, e.g. [6,7] (for an excellent review see [5] ), have not been used, in order to preserve full flexibility in the simulation of the sophisticated mechanisms of maturation in the bone marrow and premature death of newly formed red cells. Automatic choice of step length in the integration to obtain a desired accuracy has not been implemented in the programs because of the relative ease to make this choice manually in this kind of problem. In order to compare the results obtained by the methods of numerical integration with analytical solutions, the simple but "incompatible" model of Garby et al. [ 1] has been used. In the present paper this model will be called model I (see fig. 1). The "compatible" model IIc according to Vuille [3] will be called model 2 (see fig. 2). The interpretation of the symbols is given below. It is to be noticed that the exchange between the red cell compartment C and the S-compartment is very slow (the normal life span of red cells is c. 120 days) and has not been considered in the model that concerns the evaluation of rate constants. For further details, the reader is referred to the previous papers [1,3,4]. The process of premature death of newly formed red cells deserves some further comment. This process occurs both in normal individuals and haematological disorders [8-10], but the detailed mechanisms leading to early destruction of cells are not well known. It seems reasonable to assume that the red cell precursors destined for premature death represent a population that differs in certain respects from "normal" cells already at the

~ ' ~ ~ 5 ~[~4

FigureI

I---'l

L~

F-i-I

,I I-_ I, ' - - I

.___.

b

I

~

r= s

l,.----IVlh ~ --

I'.I,

L,.i Erythrocyte life span '

,'

ModelI

.o,., ~

-d

}

3 /

L___-I

tO

-o

a

Fi~lure3 Mode[ 3 Figs,I-3. Models1-3.Forexplanationof symbolsseetext, earliest stages of maturation. If this assumption is correct, this population of precursors should be separated in the model from the population of cells that survive normally. Such a model, including two different populations of red cell precursors, will be treated as well in the present paper (model 3, fig. 3). The following variables will be used: : Time after injection of the tracer. : Time variable, indicating the time of delivery to the circulation of a certain set of precursors. E(t) : Amount of radioiron excreted from the body and/or fixed irreversibly in storage. S(t) : Amount of radioiron stored reversibly in liver, spleen etc. P(t) : Amount of radioiron in plasma. C(t) : Amount of radioiron in circulating red cells.

t u

92

T. GROTH et al.

lr(t,v)

: Amount of haem-radioiron (at time t) in the precursors that will be delivered to the circulation (C) at time v. #(t,v) : Amount of exchangeable non-haem-radioiron (at time t) in the precursors that will be delivered to the circulation (C) at time v. ~(t,o) : Corresponding quantities for precursors p(t,v) that will die prematurely and deliver their iron content directly to the plasma. (In the following called "sick" precursors in contrast to "normal" precursors.) Q " Maturation time for normal precursors. : Maturation time for "sick" precursors. 60(t- v) Distribution function for the iron flow from P to p. Note that t - v is the age of the corpuscle, being - Q when the precursor is "created" and 0 when the precursor enters the blood stream. It follows that 60 ¢ 0 for arguments between - Q and 0. ff.(t - v) : Corresponding distribution function for flow from P to/a. 60 =/=0 for arguments between - Q and O. R= : denotes the fraction of the iron flowing f O ~ ( t ) d t from P to the sick precursors. Moreover, -Q it follows that f?5 60(t)dt = 1 - R . Ki : Rate constants. Subscripts according to the numbers of the arrows in fig. 2 and fig. 3. K 5 and K 7 are in the present context assumed to be identical for "normal" and "sick" precursors, though this is not a requirement for the theoretical analysis. In the simplified model 2, no distinction is made between # and if, or between rt and ~'. Instead, a fraction R of the precursors ready to leave the bone marrow are destroyed, and their iron is returned to P. In this case, we obtain the following system of differential e.~quations" dE dt = K1S"

(2.1)

dS d-7 = K 2 P - (K1 + K 3 ) S "

(2.2)

dzr d--t = K7#"

(2.3)

d/.t _ K4 p 60(t- v)- (K 5 + r7)/a • dt

(2.5)

dC dt - (1 - R ) [#(t,t) + 7r(t,t)] . dP dt

(2.7)

- K 3 S " ( g 2 + K4)P + R [p(t,t) + 7r(t,t)] + ,t + Q +/(5 Jt #(t,v)dv.

(2.8)

Model 3 yields the following differential equations: dE d--t- = g l S "

(3.1)

dS d-t- = K2P" (K1 + K 3 ) S "

(3.2)

dTr = K 7 P . d~- = K7/~ d--'t '~"

(3.3) (3.4)

dtt d--7

_

- K4 P 60(t- o) - (K 5 + KT) # .

(3.5)

(K5 +K7) ~.

(3.6)

dff d--[ = K4 P C ~ ( t - v ) -

dC d---t = #(t,t) + n(t,t) .

(3.7)

dP = K3S-(K 2 +K4)P+ff(t,t)+~(t,t)+ dt + K 5 [ f ; +Q #(t,o)dv + f;+Q-ff(t,o)do] .

(3.8)

It should be noted that if 60(t - v) = aCo (t - o), where a is a constant, the systems (2) and (3) are equivalent. In this case, we must have Q = Q. The systems of differential equations (2.1)-(2.8) and (3.1)-(3.8) cannot be solved to yield analytical expressions for the amount of radioiron contained in each compartment at any time t(as was the case for model 1, of. [1]). The systems may be solved, however, by conventional methods of numerical integration [5], e.g. Runge Kutta algorithms of different orders. Suppose that for a given tn we know the values of E, S, C and P, and that the values of/a and Ir for this t and all o ---p k ( p is an integer) are known. Moreover, suppose 6o can be evaluated for any argument (according to an arbitrary function specified as input data). Checking through the systems of equations, we easily verify that the derivatives with respect to t of all involved quantities can be computed. Using method B, the handling of the parameter v is

COMPUTER SIMULATION OF FERROKINETIC MODELS

///

t

Q

V

Fig. 4. The functions t~, ~r, f f a n d $ are defined to satisfy the conditions #(0,o) = n(O,o) = g(0,o) = ~(0,o) = 0 and #(t,t+Q) = ~r(t,t+O) = #(t,t+O.) = ~-"(t,t+O.) = 0 (fat lines).

t

tn 0

K 2K 3K

Q

Fig. 5. To get an optimal compromise between accuracy and computer time, the intervals k are made greater than the intervals h in the t-direction. The dot represents the point (nh,nh), which is found by interpolation.

I

KL 2K 3K ,K tnl

6K ,K I

tn.Q

Fig. 6. The shadowed area represents the integral Q ~(t,o)do . t The heavy dots indicate k n o w n values of/~.

a

b

c d.a b

cd÷a b

c

d

Fig. 7. The weights used in the fourth-order quadrature formula for k = Q/9.

93

tuation; the case #(t, t +Q) the initiation of new precursors. A horizontal cross section in the figure describes the state of the # pool at a fixed time t, for varying time of delivery o. A vertical cross section, on the other hand, describes the development through t of precursors with a certain time of delivery o. Referring to the system (2), we see that eq. (2.5) describes the development of # along such a vertical cross section. Eqs. (2.7) and (2.8) utilize the value o f # on the line t = u. Eq. (2.8) also uses the integral o f # over a horizontal cross section. In the procedure of numerical integration, # has to be handled in discrete points. In order to obtain an optimal compromise between accuracy and computer time it is desirable that the intervals k in the o direction are greater than the intervals h in the t direction. This condition is shown schematically in fig. 5. At a fixed point tn = nh we may know the value of # in points (k, nh), (2k, nh), (3k, nh), .... but due to the difference in interval length, we do not know explicitly the values of # in the point (nh,nh) which is used in eqs. (2.7) and (2.8). Of course, we can obtain the value of #(nh,nh) by interpolation in the v direction. In order to make this interpolation efficient, we have to integrate # by (2.5) beyond the line t = o. The vertical lines in figure k indicate how far we integrate. It should be noted that the value of #(t,o) for t > o has no physical significance, but that it is perfectly possible to use it in the computations. In equation (2.8) we need

f

t+Q #(t,v)dv.

t Considering a horizontal cross-section for t = t n in fig. 5 we obtain the desired integral as the shadowed area in fig. 6, where heavy dots indicate known values of #. A special fourth-order quadrature formula has been designed for the calculation of the desired integral. Let A be a measure of the skewness, defined by t = ( p + A ) k , where p is an integer and 0 ~< A < i. For example, A = 0.5 for tn = 3/2k. Next, define a = 3 + 6 A 2,/3= 8 A + 4 A 3 , 7 = 9 - 6 A 2 , 5 = 12A3,

not quite trivial. The functions # and 7r are defined on the strip 0 ~< u - t <~ Q (see fig. 4). Restricting the discussion to #, we have #(0,o) = 0 and #(t,t+Q) = 0 (fat lines in fig. 4). The case #(0,o) describes the initial si-

(3.9)

94

T.GROTH et al. Table 1 Parameters used in the computations

and a=

(3.10)

Parameters kl

it is easily verified that the formula v+3h+Ah

ka f ( t ) d t ~ af(v) + bf(v+h) + cf(v+ 2h) + ks k7

+ df(v + 3h)

Model 1

Model 3

0.0000 0.5482 0.0555 9.2478 0.2064 0.2920

0.05 1.50 0.15 7.50 2.00 9.00 0.100 7.50

R

(3.1 l)

gives the exact result if f is a polynomial of order ~< 3. By iterated use of this formula we can obtain good approximation of

a Q

6

ta)

Rectangular distribution function

Triangular distribution function

t+Q / t

#(t,v)dv

if Q/k (the number of steps in the v direction) is a multiple o f 3. For instance, if k = Q/9 we must use the w~,'ights indicated in fig. 7. (Note that d + a = ¼a. The function 7r and the corresponding functions ~ and are handled in the same way as #. As a check of the integration, the total sum o f all the compartments is used ir~ every timestep (conservation of amount of injecl:ed radioiron). Deviations greater than a specified size are printed out as error messages.

2. Basic flowcharts The basic flowchart for method A is shown in fig.8. This flowchart applies to all three models in general. The bone marrow processers are, however, illustrated for model 1 and the last box in fig. 8 has to be modified in order to cover the details of models 2 and 3 (see figs. 1,2 and 3). As a detailed complement and in order to illustrate the simplicity of method A the corresponding F O R T R A N program for model 3 is listed at the end of this paper. In fig. 9 a brief flowchart is given for method B as applied to model 3.

3. Typical sample runs Sample runs will be given at the end of this paper for methods A and B. Model 3 has been chosen. Table 1 gives the parameters used in the computations.

Q Q1

Q2

01 = 0.5 × Q2 02

f

to(u)du = 1

0

The theoretical radioiron curves in compartments P (disappearance curve) and C (appearance curve) are shown in figs. 10 and 11. As a visual illustration of the rate of convergence the disappearance curve is given as calculated with methods A and B using various step lengths. From these figures it can be seen that: 1. Both methods converge towards the same limit curve.

2. The Runge - Kutta method converges, as expected, more rapidly than the simulation method. In fact, the differences between the curves obtained with the Runge - Kutta method and step lengths from 0.001 to 0.05 day were so small that they could not be drawn in the figure. 3. The deviation between the two methods and between two curves with different steplengths is largest within the interval 0 - I day, while the differences are negligible after t = 2 days. The problem of accuracy in the integration can of course be managed by automatic step length control, but in order to reduce computer time needed in the calculation of a large num-

COMPUTER SIMULATION OF FERROKINETIC MODELS

I READ INPUTDATA, NUMBEROF BONEMARROWSUBCOMPARTMENTS TIME STEP DT RATE CONSTANTSKij MATURATIONTIME O. TOTALLENGTH OF CURVES TOBE PRODUCED,TMAX

I [EAcVAoLUA:EUTPTAKEE RIATTEI;OuRI~E

NT

SET INITIAL CONDITIONSI ALL POOLS[XI)=0 TIME VARIABLE T=0

i

I

'INJECT'"RADIOIRON,(TOTALAMOUNT,ill INTO PLASMAPOOLP

"I INCREASE T BY DT I

I CALCULATE AMOUNTOFRADIOIRON FLOWING FROM COMPARTMI TO COMPARTMd DURINGTHE INTERVALDT FjI" KjI' XI • DT

i

IJ

ADD AMOUNTSFLOWINGTO X] SUBTRACTAMOUNTSLEAVING X!

i DISTRIBUTEAMOUNTFLOWINGTO PR / AMONGSUBCOMPARTMENTSACCORDINGI TO THE RESPECTIVEUPTAKERATES /

I

LET SUBCOMPARTMENTSPROGRESS I STEPALONGTHE Q-AXIS L,- I .

I.

•~ I • 1 i

I/

I/ ~

I. "

I,"

I

] • ~?:ii::i::ii!]

ADO CONTENT I OF LASTSUBCOMPARTMENT CREATE TO CIRCULATINGBLOOD NEW SUBGOMP Fig. 8. The general principles of method A presented as a block-diagram.

95

96

T.GROTH et al.

READ INPUT DATA AND DEFINE:

RATE CONSTANTS TOTAL LENGTH OF CURVES TO BE PRODUCED: TMAX MATURATION TIME OF PRECURSORS: Q TIME SPACING OF SUBCOMPARTMENTS IN THE BONE MARROWs DV TIME STEPs H (DV-Q/NI! H=DV/N2! NI, /12 INTEGERS AND NI-MULTIPLE OF 3)

SET INITIAL CONDITIONS E_-S~C-O, P-I ALL S U B C O M P ~ N T S ~ -x.U=n:O TIME T-O

("INJECTION" OF RADIOIRON)

I FIRST R-K STEP: EVALUATE AND

t-T t+Q / ~(t,v)dv AND t

I t

~ ( t , v ) ~ ACCORDINGTO (3.o9-3.11)

u(t,t), x(t,t)j~(t,t), ~(t,t)

BY INTERPOLATION

EVALUATE THE DERIVATIVES dE

dS

dP

dC

d,

d,

d~,

d~

ACCORDING TO (2. ~-2.8) EVALUATE NEW VALUES OF E, S, P, C, ~, ~ AND x, ~ FOR TIME t+H/2 BY PROCEEDING HALF A TIME STEP ALONG THE TANGENT. SECOND R-K STEP,

t-T+H/2

REPEAT EV£LUATIONS IN THE FIRST STEP

THIRD R-K STEPs

t-T+H/2

REPEAT EVALUATIONS ABOVE ,PROCEED H

FOURTH R-K STEP,

t-T+H

REPEAT EVALUATIONS ABOVE EXCEPT THE LAST POINT.

THE VALUES OF E, S, P, C, ~, ~ AND x, ~ ARE EVALUATED IN t-T+H FROM THE CORRESPONDING VALUES IN t-T AND A WEIGhtED SUM OF THE DERIVATIVES IN THE FOUR R-K STEPS. E.G.

sCT+H) - S(T) + HiS

4Z , % ( ~ ) i i=I

wi'I,2,2,1 T-T+H

PROCEED ONE TIME STEP

YES

L ~ S~OO_~,aTMEnS I ~, ~, ~, p PROGRESS[ I STEP DV ALONG THE I Q-ArtS. CHEATE NEW I SUBC OMPARTMENTS. J

NO NO

Fig. 9. The general principles o f m e t h o d B presented as a block-diagram.

I

COMPUTER SIMULATIONOF FERROKINETIC MODELS % radioiron in plasrna

Curve Method Step number length 1 Runge-Kutto 01250 2 - "0.05-OD01 3 Simulation 0.g025 4 -"0.~ 5 -"0.0125

100.

~I 50 i 40. ', 30 ,,

5

-'-

oo~

7

- - -

~0s00

8

-"-

o.tooo

97

% radiairan in p l a s m a 1.0-

20'

1 0.5- 2 -

10.

0.35. 4.

0.2-

3. 2.

0.1-

1-

x .---,, '~0.1 ;"-'-36

0.5-

Og5-

0

days after injection

I

2

3

4

5 6 7 8 9 10 11 12 13 14 15 16 days offer injectian

Fig. 10a. The disappearance curve as calculated with various step-lengths and parameters according to table 1. t = 0 - 0.8 days. Fig. 10b.Ttie disappearance curve as calculated with various step-lengths and parameters according to table 1. t = 1 - 16 days. Curve numbers as in fig. 10a. %

Table 2 Accuracy of the two methods on the basis of model 1. The disappearance curve. "True" curve obtained from analytical expression.

radioiron in r e d c e l l s 100

0f 0-< / 0

1

.

2

.

.

3

.

4

.

.

.

.

.

.

Method

Step length (days)

A

0.005 0.002 0.001

5.5 3.0 1.5

0.01 0.05

4.5 0.2

.

5 6 7 8 9 I0 II 12 13 14 15 16 d a y s after injection

Fig. 11. The appearance curve.

ber of similar curves (parameter estimation) it is suitable to choose two different step lengths and switch the step length after a certain time. In the present case, t = 1 day was found to be optimal for switching the steplength. The optimal "switch-point" may be different for different models and different parameter combinations, but can easily be found empirically. The convergence properties of the two methods can suitably be investigated using model 1, for which "true"

B

Maximum deviation (%) from "true" disappearance curve.

reference curves can be obtained from the analytical expressions in reference [1]. It was found that with decreasing dt both methods yielded curves that in every part converged towards the "true" curves. For the disappearance curve, the maximum deviation from the true curve is shown in table 2. Table 3 shows the corresponding results for the appearance curve; for this curve, the accuracy is defined as the deviation at t = 1 day, since the maximum deviation always oc-

98

T.GROTH et al. Table 3 Corresponding data to those in table 2 for the appearance curve Method

Step length

A

0.005 0.002 0.001

2.0 1.5 1.0

B

0.1

0.2

Deviation (%) from "true" appearance curve at t = 1 day

Data Center, Uppsala, Sweden. The computer is a 32K, 48 bits/word computer equipped with 2 magnetic drum storages (524K each) and 6 tape drives. Storage requirements and time for sample runs: Model 3: Method A ~ 10K, c. 4 sec. Method B ~ 5K, c. 90 sec. The storage requirements are highly dependent on the number o f subcompartments in the bone marrow compartment.

5. Mode of availability curred at t < 1 day, a region o f the curve with no practical importance. From these tables it can also be seen that method B converges much more rapidly towards the " t r u e " curve than method A, requiring considerably less computer time. With the conditions that 2 = Q and ~ = co, model 2 and model 3 yield identical curves. Because of the double number o f subcompartments/a and n in model 3, however, the computer time needed for the production of one s e t of curves is increased by a factor o f almost 2. Therefore model 3 should only be applied in cases, where the experimentator has got a wellfounded hypothesis concerning the kinetics o f maturation o f "sick" precursors. In this case, the same conclusions as regards computer e c o n o m y o f the methods A and B will apply, as were obtained on the basis o f model 2.

4. Hardware and software specifications All the programs are coded in F O R T R A N IV for the CD 3600 computer at the Uppsala University

The programs are available from Uppsala University Data Center Box 2103, 750 02 Uppsala, Sweden.

References [1] L.Garby, W.Schneider, O.Sundquist and J.-C.Vuille, Acta Physiol. Scand. 59, Suppl. 216 (1963). [2] MoPollycove and R,Mortimer, J. Clin. Invest. 40 (1961) 753. [3] J.-C.Vuille, Acta Physiol. Scand. 65, Suppl. 253 (1965) 3. [4] J.-C.Vuille, Acta Physiol. Scand. 65, Suppl. 253 (1965) 63. [5] W.Jentsch, Digitale simulation kontinuerlicher Systeme (R. Oldenbourg Verslag, Mtlnchen, Wien, 1969). [6] M.Berman, E.Shahn and M.F.Weiss, Biophys. J. 2 (1962) 275. [7] Y.Chu, Digital simulation of continuous systems (McGraw-Hill, New York, 1969). [8] L.G.Israels, J.Skandenberg, H.Guyda, W.Zingg and A. Zipursky, Brit. J. Haemat. 9 (1963) 50. [9] I.M.London, R.West, D.Shemin and D.Rittenberg, J. Biol. Chem. 184 (1950) 351. [ 10] T.Yamamoto, J.Skandenberg, A.Zipursky and L.G. Israels, J. Clin. Invest. 44 (1965) 31.

COMPUTER SIMULATION OF FERROKINETIC MODELS

Program

listing

Method

99

A:

PROGRAM MOOrL 3 IRON-FLOW IN TH~ HUMAN BODY,MFTHOD A TMAX=TOTAL LENGTH OF CURVES TO BE PRODUCED READ INPUT DATA TSHIFT=TIME FOR CHANGE OF TIMESTEP DTI~TIME STEP BEFORE T=TSHIFT DT3=TIME STEP AFTER T=TSHIFT NF=NUMBER OF SUBCOMPARTMENTS IN BONE MARROW POOL DT2=TIME INTERVAL IN WHICH B.MARROW POOL IS DIVIDED **OBS** DT2=P2/NF e*OBSe* DT! AND DT3 LESS OR EQUAL TO DT2 DT21DT1 AND DT2/DT3 MUST BE INTEGERS P2=Q IN PARAM,LIST=MATURATION TIME FOR PRECURSORS Pl DFFINFS THE POSITION OF THE TOP OF TRIANGULAR DISTRIBUTION FUNCTIONS **OBS** INTHIS PROGRAM PI,P2 ARE THE SAME FOR SICK AND NORMAL PRECURSORS.WHEN OF INTEREST THIS MAY EASILY BE CHANGFD S I , , . . , . . , S T = R A T E CONSTANTS S6=PART OF IRON FROM PLASMA TO COMPARTM. SICK CELt. I X , I Y ARE PRINT INDICATORS, PRINTING OCCURS EVERY IX THIBEFORE TSHIFT) RESP, IY TH(AFTER TSHIFT) STEP DIMFNSION AFACH(]OO0),DD(IOOO),EE(IO00),AFACH2(IOOO),DD2(IOOO),

C C C C C C C C C C c c rtfC C C C

1~2(1000)

60 P ~ ( 6 0 ) 1 0 0 2 } SI)52,$3)S4)$5)$6,$7 W~IT~(K,I,1020) 1020

FnPMAT(1H1)

WRITF(61tIOO~)SItS2,S3,S4,SStS6,57 ]004 FOP~AT(SX,*RATF CONSTANTS*/(3F20.IO)) IFIST)61,61,6P 1OO2 ~ORMATIBF]O,O) 62 PF#O(60,lOOl) I X , I Y WRIT~(&],IO05) I X , I Y 1005 FORMAT(SX,*PRINT INDICATORS*t2110) RFAb(60,IOOO)TMAXtDT],DT2,ISHIFTtDT3 Inno ~RMAT(BFIh.n) WRITF(61,1007)TMAXtDT19DT2,TSHIFTtDT3 1007 F~RMATISXt~TIME PARAMFTFRS~/(3E20.]O)) PCA~(6Og]nOR)NC,PIgP2 ]nh8 F~RMAT(IIh~2FIn,O) WRITF(6I~IO~q) NF,PI,P2 10~9 FORMAT(~X,*DATA FOR DISTR,FUNCTION*gIIOt2FI0.2) WRITF(61tZO11} 1011F~RMAT(////) WRIT~(61,!OO6) ]006 FORMAT(SX,*IRON-FLOW IN THE HUMAN BODYtDIRECT SIMULATION,MODEL 3*) 1OO1 ~RMAT(3110) WRTT~(61,IO]~) 1010 FORMAT(6X~3H T ,&XtIHE~TX,SH S ,4X~6HPLASMA,~X~8HB.MARROW~2X, ]SH~I.~hD ) nF~tN~ THE ~ I S T R I R U T r O N B~ 12 I=I,NF

AI=I I~(O-~I)10,10,11

FUNCTIONS,HFRE

OF

TRIANGULAR

FORM

100

T . G R O T H et al.

10

AF^CH(TI:2.~P~DT2/(P]*P2}*(I.-S6) AFhCH2(II:2.~P*~T2/(PI*P2)*S6 ~n T~ 12

II 12 C m N T T N U F

~VFTN~ INITIAL CONOITTONS OF ALL POOLS A:POOL E IN FIGuRV R:OO~L S C:PO~L P D=POOL M, DD=SUBCOMPARTMENTS IN D F=POOL PR, EE=SUBCOMPARTMENTS OF E DE=D+E= THE BONF MARROW COMPARTMENT INDEX 2 FOR SICK CELLS,E2~D2,EE2,DD2

FG=POOL C DV=INF TIMV,TIPESTFP,TIME FOR SHIFT OF SUBCOMPARTMENTS IN B,MARROW POOL TNJFCT RADTOTRON (TOTAL AMOUNT=If INTO PLASMA POOL ~n 2n T:I,NF DD(T)=O. D~2(I)=h, VV2(I):h. )n FC(T)=Q. ~T=~TI

IVFR=~T2

C:I.

r):~. F:~,o

D2m'~o F:O,

~V:O,

DF~:O. PhQLS:I. NP=~ AUS=I.-(SS+ST}*~T ~IN=STe~T

~ IF(NP-NRITX~TX)~O,32,40 ~2 WRITF(61,10n3)T,A,B,C,BENM,FG 1003 FnP~AT(FIno2,6F]O,6| I~(ABS(POOLS-I.)-I.E-6) ~0~40t42 a2 WRTTFI61,1~I2)POOLS Ih12 FORMAT{SX,~OBS SUM OF ALL COMPARTMENTS=*,E20.10) INCREASE TIME BY DT ~

T:T÷~T

RFIN=C*S4~DT CAICULATE AMOUNT FLOWING FROM AND TO THE VARIOUS COMPARTMENTS AND ADD -SUBTRACT FROM THE POOLS DISTRIBUTE AMOUNT FLOWING TO PR AMONG SUBCOMPARTMENTS ACCORDING TO THE RFSPFrTIVF UPTAKE RATES D~

al

I=I,NF

FV2rI)~EE2fII+O~2(I)*EIN DD2(1)=DD2(II~AUS+AFACH2(I}*REIN FF(II=EE(II+OD(I}*~IN

COMPUTER SIMULATION OF FERROKINETIC MODELS

10

AI DD{I}=DD(1)~AUS+AFACH(1)~RC!N F=F+D*~IN F2=F2+D2*EIN DALT=D DAIT2=D2 D=D~AUS+RFIN~{1,-S6) D2=~2*AUS÷RFIN~S6 A=A+B~SI~DT RALT=B B:B+(C*S2-R~(SI+S3))*DI C=C+(DALT~SS+DALT2*SS+BALT~S3-C~(S2+S4)I~DT 4 4 IFfTVFR-T)A6,4&~b2 LEI SUBCOMPARTMENTS PROGRESS 1 STEP ALONG THE O - A X I S . ADD CONTENT OF LAST SUbCOMPARTMENTS EEtDDoEE2,DD2 TO CIRCULATING BLOOD AND PLASMA=CREATE NEW SUBCOMPARTMFNTS 4 6 TVER=TVER+DT2 470 FFALT=EE(I} DDALT=DD(1) FFAt.T2=EE2(]) D~ALT2=DD~{I} 47 D~ 51 I=2,NF Fr(I-I)=EE{I) FF2(I-I)=FF2(1) DD2(I-I)=DD2{I) 5! DD(I-I)=DD(1) EE(NF}=O, DD(NF)=O= FF?(NF)=O= D~2(NF]=O, ~ D=D-D~ALT D?=D2-DDALT2 F?=F2-EEALT2 F=E-FEALT F=F÷DhALT G=G+RFALT C=C+FEALT2+DDALT2 ~2 DE=D÷F DF2=D~÷E2 BFNM=DE+DE2 FG=F+G POOLS=A÷B+C+DF+FG÷DF2 N~=N~+I TIMESTEP DT IS CHANGFD WHEN T EQUAL TSHIFT IF(ABS(T-TSHIFT)-IoE-8) 4 5 t 4 5 , S g 4~ ~T=DT~ AUS=I.,-(SS÷ST)tDT ;IN=ST*~T IX=IY NR-O 59 IF(T-TMAX)30,30,60 61 CALL FXIT FND

102

T.GROTH et al.

Method A. RATE CONSTANTS 5,0000000000"002 ~,5000000000"000 7,5000000000.000 2,0000000000*000 9,0000000000~000 PRINT INDICATORS 4 20 TIME PARAMETERS 2.0000000000"001 2.5000000000-002 ~,0000000000-000 5,0000000000=002 DATA FOR ~ISTR,rUNCTION 75

1,5000000000-001 1.0000000000-001

1,0000000000-00! 3,75

7,50

|RON=F~OW IN TWE HUVAk BOEY~CIRECT SIMULATIOk,MODEL 3 T E S PLASMA B,MARROW BLOOD 0,00 0,000000 0,000000 1,000000 0.000000 0,000000 0,10 0,000241 0~106748 0,393824 0.49902E 0,000160 0,~0 0,000868 0,149269 0,183945 0.665252 0.000665 0,30 0,00165~ 0,167~21 0,093547 0,735384 0,001593 0,40 0,002512 0,175607 0.049869 0,769026 0,002986 O,bO 0,003396 0.17~122 0,027754 0,785861 0,004868 0,60 0,004~87 0,177976 0,016105 0=796764 0=004868 0,70 0 , 0 0 5 1 7 5 0,176480 0,010189 0,800903 0.007254 0,80 0,006053 0,174316 0,007109 0.802369 0,010153 0,90 0,006920 0 17~834 0,005504 0,~02171 0,013571 1,00 0,007774 0 169~15 0,004668 0.800833 0,017510 2,00 0,015621 0 144547 0,003858 0,750359 0,086015 O0 0,022327 0 123901 0,004043 0,641106 0.208623 4 O0 0,028126 0 107943 0,004299 U=473140 01386493 5 O0 0,033204 0 094772 0e004751 0,287806 0,579466 6 O0 0,037651 0 082019 0,003646 0=169486 0.706598 7 O0 0,041513 0 071408 0,002530 0,099102 0,785447 8 O0 0,04483~ (1 060993 0,00~660 0.074176 0,818339 9,00 0,047659 n 051~77 0,001420 04059000 0.840044 10,00 0,050064 0.044128 0,001201 0,045989 0,858617 11,00 0,052109 0.037521 0,001001 0,035865 0.873504 12,00 0,053848 0,031870 0.000829 0=028515 0=884930 13,00 0,055324 0.027065 0.000690 0,023218 0,893703 14100 0,05~677 0,022966 0,0005?9 0.019240 0,900638 15,00 0.057641 0,019483 0t000488 0,016076 0,906312 16,00 0,0~8543 0,016527 0,0004~3 0,013461 0,911056 17,00 0,05930d 0,014018 0,000350 0,011295 0,915029 18t00 0,059957 0,011887 0,000274 0,009828 0.q18054 19,00 0,060507 0.010082 0,000232 0.008269 0,920910 20,00 0,060973 0,008~51 0,000197 0,006956 0.923322

COMPUTER SIMULATIONOF FERROKINETIC MODELS

103

Method B RATE CONSTANTS 5,0000000000=002 1.5000001~000-000 1,50000000(I0-001 7,5000000000~000 2,0000000000*000 1,0000000000-00! 9,0000000000-000 TI~b VARIABLES 300 i 6n 4 QmO~eT)EL, 7,5000000000~000 7,5000000000~000 1.0000000000*000 T~AX,TSHIFT 2.00000flO000*O01 1,O0000qO000*O00

IRON.FLOW I~ THE HUPA~. BOCY,MOCEL 3,RUNGE-KUTTA METHOD DV~H 2,5000000000-bn2 2,5000000000-(102 T E S ~LOCD B,MARROW PLASMA 0 000000 0 000000 0.436140 0 ~00000 0 0n0285 0,099665 U,OO00aO 0 .100000 0 463820 0,212655 0.200000 0 , 0 0 0 9 0 7 0.14~51}2 0 , 0 0 0 4 5 0 0 200000 0,642485 0,111296 0 300000 0,001681 0,163778 0,001211 0 300000 0 722030 0 060985 0 AO0000 0,002525 0,172902 0,002416 0 400000 0 76116~ 0 034723 0 500000 0,003401 0.176374 0,004096 0 500000 0 781402 0 020669 0 bOO000 0 004~b 0,176885 0,006270 0 600000 0 791888 0 013061 0 700000 {! 005167 0 175828 0,008949 0 700000 0 796991 0 008925 0 800000 0 006042 0 173047 0 012141 OobO0000 0 798942 0 006674 0o900000 0 006906 0 17~643 0 015852 0o900000 0 798921 0 005453 lo000000 0 007758 0 169135 U 020084 1,000000 0 797594 0 004221 2,000000 0,01557~ 0.144336 0 0~1568 2o000000 0 744294 0 004605 3,000000 0,022272 0 124%65 0 217270 3,000000 0 63~673 0 005005 A,O00000 0 028006 0 108246 0 397485 4,000000 0,461176 0 004223 5,000000 0 033140 0,094936 0,572745 5o000000 0,294955 0 003300 6,000000 0 037~79 0 082820 0 702133 6o000000 0,174175 0,002365 7~000000 0 041437 0 071639 0 783328 7o000000 0,10%200 0,001592 8,000000 0 044756 0 06~246 0 817858 8o000000 0,074537 0,001362 9~000000 0 t]47584 0 05~130 0,839723 9o000000 0,059.200 0,00114b 10o000000 0 049992 0 044~71 0,858372 lOoO00000 0,046107 0,000957 11o000000 0,052040 0 037750 0,873294

104

T.GROTH et al.

11 000000 12 000000 12 000000 13 0 0 0 0 0 0 13 000000 14 000000 14 000000 15 000000 15 000000 16 000000 16 000000 17 000000 17 000000 18 000000 18, 000000 19, O000OO 19 ooooon 20 000000 20 000000

0,053183 0,055263

0.032092 0.027263

0,884718 0,893485

0,035945 0,000796 0,028597 0,000065 4

0,056621 0,05758~

0.023148 0.019649

0.900411 0,906076

0,0233~0 0 000559 0,019347 0 000473 0.016199

0,0~8494

0.0J6678

0,910812

0 000400 0,013601

0.059264

0.014155

0,9!4779

0 000339 0,011449

0m059916

0,012913

0t9~8~,16

0 000287 0.009663

0.060470 U.060940

0,010194 0.008650

0,920905 U.923267

0,000243 0.00B172 0 000206 0,006920