- Email: [email protected]
Finite element methods in cell dynamics fundamentals
1. lnrroduction The matD~matical analysis of many biophyrical problems often involves the solution of on,: or more differential equations Criphonov and Samankic, 1963). It may be diffint!t or impossible to obtain closed fan so!utiotts by conventionii means, and many techniques or methods have been developed to circumvent the mathematical dilficuities. Finite element analysis is at? approach which tnchider the following prmciplcs: 1) The fundamental physical or chemical lws governing B given system are used in as direct a form as possible with a minimum of mathenstical analysis or use of analytic solutiaa. 2) The variablesin the laws are handled in incremental forr) and are obtamed by summing over the vxiable incremwts as defined by the laws.
There approaches are in contrast to tbc classical methods of (a) mathematical analysis and(b) convcntional computer programming. Commonly, the laws governing a system are analyzed mlthematically by dcrtving differential equations and by solving them in closed form, wherever possible. These differential equations, or particular solutions of them, are ttad at the beginnir.g 0; suctt a mathematical analysts. In contrast, finite
R
MOORE. N.DAYIDS
and RLBERCER
studied. This may be accomplished by using 8 running index t” designate which parts oft problem arc being considered. The physical law, in incremcnta’l form, is tberl staed. Vsrious examples arc pen further in the text. Tbc initialization of the variahlrs is necessary t” rtiot tbc cycle ofcom~utations, and a set of bandary conditions is needed which arc similar to those in clas:ical mathematical analysis. In the example of a freely falling body, to be discussed. if a weir&t has ari initial velocitv of zero. and if it falls from a heaht of one hundred feet, the”: u= O,andy = Ior,, for r=O. An important part of the calculations ia fimte element analysis involves incrementing the variables appearing in the physical laws and summing over these increments by looping back in the program. The magnitudes of the increments arc cb6s:n to satisfy criteria of soccitied accuracy and “erhaps #criteria of stability. In a-given problcm,wrt& ratios, referred to a~ “moduli”, arc useful in f&ing the ratio “f the incremcots, such as that between distance end time. The next step in the constmctit~n of a finite ~clcment analysis program is to establish the loops by which the summation over the vaoables is performed. A “loop” consists of ““c or m”rc
xe
.
F&z. i. scieme for problem soI”¶iO” by ahnate
methodi.
In conventional approaches to the solution oim:tthcmatical forms of biophysical problems, the invcstigator may start from a previously derived mathcmatic.d solution in closed form. expressed a an mtegr;d “I as a series. Hc may ultimately apply a computer to helo him. A “finite difCeerence” method. or other nunwical method, is frequently used to arnve at a solution. w!wn such a closed-form solution is unavailable. Finite clement analysis, by cootmst, operates on tbc laws themselves, not on derived relations. The following steps arc necessary in carrying out the finite element attalysis of a pm5lem: i) statcmcnt of the probLcm, 2) subdivision of the physical structure into finite elcmcna, 3) initialization of the variables, 4) calculation in sequence (by “loops”), a.. exprcssirl of the physical laws in incremental form. b. the summa!ion of the wriablx The statemen: of the problem mus? be phrased in ss specific terms as possible. Thus, the transport process thraugh cellular membranes may bc resolved into a :argc number of specific components. It is useful to desigr.ate the specific components which arc being
3. Averagingof initial and final velocity over the time interval: ” = ;(ui + q). 4. (a) Incrementing the height: Ayv=uAt. (b) Accumulation of height mcrementr, yf =yi + Ay. 5. We then mcrement tho time, 1, by a step. AI, and repeat the above cycle of computations. In this simple case, the force is assumedto be constant. If the force is not constant, as in simple harmonic motion, then correction terms have to be used to calculate the final velocity. The discussionof this problem occurs ekewhere (Balko and Berger, 1969). A suitable magnitude for the time increment can be establishedby trial. The problem is mn with EUCcessivelysmaller time increments, until the graphs of successivesolutions superimpose.Once the solution is found to be invariant with respect to choice of time increment, the value for the time increment may bs mcreasedto conservecomputer time, but still to obtain a solution within the desired bounds of accuracy. Thts pant will be dealt with in depth in later papers. The solution to a simple problem is: A” =gAf = 32.2 X 0.1 = -3.22 ft/sec , uf = u, + Au = 0 - 3.22 = -3.22 It/xc u=t(ut+u,)=:(O-3.22)=-I.61 2. Simple application
Ay==uAr=-1.61
X0.1 =-0.161
, ftlsec,
ft,
yf =yi + Ay= IO0 - 0.161 = 99.R36 ft. This section will describe the application of the finite element approach to a simple problem, the mu tion of a freely-falling body. Later (Moore, Gwids and Berger, 1969; Balk0 and 13erger,1969), more complex problems, of interest to biology, will be dedt with, since the greatest advantageof the method is its apphability to highly complex systems. 2.1. Andysis of mor;oon in c gmvimtioml fidd :. Initialization ofvariables: At the beginnmg softhe problem, the time, I = 0, and the velocti~i. ” = 0, and the = IOQ feet. 2. Incrementing the variables: (a) the velocity iwement, aV =gAt, (b) the accumulation of velocity is done by addmg the increment, Au, to the initial value ~oCvelocity. u,, to get the final value of velocity. q: uf=ui+Av.
height,y
rf=ti+Ar=O+O.l=O.lsec. Ii the redder cares to perform this set of calculations for twenty-five times, he will obtain a table of results similar to the output of the computer program. 2.2. Computerpwgmnt The jldzated cumulation of the variables may be performed by machine. This is done by including these operations in a loop. A loop is a set of instrw Lionswhich are executed one or more times. The cumber of times the calculations in the loop are repeated may be controlled in severalways. This section will discussthe kinds of “BASIC” statements that may be used to control a loop: 1) “FOILNEXT” statements 2) “LF...THEN” state!ne”ts 3) Independent mop counter
R.MOORE.N.DA”IDS
98
A “FOR...NEXT”’
statement is a pair of instmc-
tions of fhc form: sI=IT03
STEP I
I i!i tbc indcpcndent variable; the “FOR” statement MS its value to I, and the subsequentinstructions are perFormed down to the “NEXT I” statement. The program goesback to the “FOR” statement. increases the xdue of the I by i from the instructions to the ‘NEXT” statement, goesback to the “FOR” slatemer.t, and increasestbc value ofI by from 2 to 3, at which point I equals the second nunbcr in the “FOR” statement. The instructions arc performed for the last time. curdthe prcaam continues prist the “NEXT I” statement, to the remainder of the pmgratn. In the pr5grun cntitlcd “GRVTY I”, “FOR...NFXT” statanents ue used in lines 520 and 6 IO. Fig. :#showsa conpanson between the analytic and the finite slemen: solutions.
I tr,2,exewtes
1
and R.L.BBRGER
condition to tc met. If the condition is met, then the prol~zm ;s directed to a line, identified by a number, to cw tinue the calculation. If the con&i!on is not met, then the next ins:ruction in !equencc :s executed. For c:xnmplc,“IF T < 2.6 THEN .jZo” will cumparc the a rrcnt value of T with 2.6. If T is less than 2.6. then the comptttcr goesback to Ii le 520 and ~‘xewtcs that tolstmction.1n”GRVTY 2”,“IF...THEN” stcbcmew! are used in lines 520 and 6’10. In “GRVTY 3”. an independent loop counter is employed, rather than looping ba :cd an the valuesof any OFthe val-iables,as in lines 520 and 590. It should he noted that the symbols Dl,DZ,... are intendsd to stand for At nnd Ax. becauseof the unavailability oi the Greek “A” symbol and the limitations in naming I variable in the “BASIC” langual;c. However, D I. D2,.. should not be confused with differentials. The sectitmsof the program arc mtmb’eredand labeled in which major portions of the work ore done. The subdivisions of the computations section, by groups of instructions, are shown in the table below.
F%yri:allaw: Initkilting tll~ “1riable-. tllCt~W”ting tie “ti*bt:l Lwp controt
The solution was checked by calculating the six _ displscemcntsfrom tie cxact formula y = 100 -OS&. f = 0.1,0.2, .. .. 2.5. The results agreed to three decimal placer. An “IF...THBN...” statement is one which states a
530.s70
53O,S?O
550.580
310-370
3to-370
310.-3,0
540,560 520.610
520.540.580 520,6,0
540.S60.590 5ZO.610
III other words, the three gravity programs illustrate due: methods of looping. In I, the looping is controlled by the instruction, “520 FORT = 0.1 TO 2.5 STEP Dl”. This setsthe variable T to an initial value of0.l. executes the instructions down to “610 NEXT T”, and returns to line 520. The value of T is increasedby the amount DI (“331) LET 01 = 0.1’3, ar-i the instructions are again executed to line 610. Thi:. is done until T= 2.5, for which the instructions arc xccuted fur the final time. In II, the control of the loop is different. and it is done by”610 IFTC 2.6THEN 52O”.Tbe initial ~~al~eofTI;(T=T+DI=O+O.l)=O.i.Theinsrructior s are executed down to line 610. This compares
the value of T wth 2.6. If the value of T is less than 2.6, the program returns to line 520. If T is equal lo 3.6, then the program stops. In 111, the looping is controlled by line “530 FOR I = 0 TO 24”. The variables HPCsubscripted, and the value: of the subscript, I. increases by 1 for rach time incren~t, until the value of the subscript equals 24 Then the instructions are executed for the linal time. I is, in effect, a loop counter which keeps tlack of the number of times the program goes thm the set of calculations in lmes 540-6GO. There are additional kinds ofcontrol of lowine possible. One may have several nested loop:;. l&g combinations of the above three kinds of control. and there may be severa: nested loop counters. There may be comparison statements for variables othsr than the independent variable. T&se may be used so as to causr the program to generate families of solutio~a. for many ranges of the variables, and even to search for, and to print out, only those few solutions which are of interest for some particular reason. The availability of several means of conrrollmg loops frees the investigator from limitations inherent in either the notation employed or the mbrcnpts allowed by the language.
3 Compartmental
where 01 = change of amount of substance in the corn. partmen,, IL= rateconstant for the process, A(1) = amount of substance in the comoartmrnt at the Itb time. D2 = time mcrement. After DI has been caiculated. the increment ofA is computed, and this increment is added to the previous value A. This new value computed for A is ~tilbed m the next increment reque&This is accomplished by the instruction “LET A(l+I)=AlIPDl”. ‘Tlme is then incremented, using the in&c& “LET T(l+ I) =T(I)+DZ”. This also performs the operation of loop ing back for successive increments. The program then calculates A(0 versw r(r) for three different values of the first order rate constant. k, as shown in fig. 4.
of
analysis, part 1
Compartmental analysis may be used to describe the movement or change of substances in luantitattve fashion (Hart. 1963). For example. finite clement analysis enables the reduction of data from radioactive tracer studies to be performed in a very simple WY.
In a chemical :eaction in whxh A gces to B, v&ten A&B, the rate equation governing this process, called first order kinetics, is given by &I = -kAAi
:.te
where bA is the change in amount, A, in cornpartment in the time At; k is the first order rate constaM inset-t. In “BASIC” notation this is writ’en as: DI = -l?A(I)*DZ,
The variables are handled in the following way: time is considered the independent variable, and the amount of substance present in a cotr.partment is the dependent variable. The amount initially present, and the rate constant, are program parameters which may be changed to fit various experimental conditions. 3.3. Dmmenrarion The documentation perts:
consists of the following three
cially one which is B function of :he input rate), and coup!ing between lass rate and input rate by mwlr of wriox functions (with and witbout time delays). Each casxcan be explored by changing just ora: to three lines of the program.
4. Compartmental analysis. part :! (prugmm name: COMP2)
:;tatemevr ofrhe
4.1. problem The problem is to determine the kinetic bebaviaur of a ‘we-compartment system which cxchaogescontents by first order Proresser(see fig. 5).
3.4. Discussion This case is simple but common. It describesrsdi* isotopic decay, first order chemical reactions, operation of a shock absorber, cooling of a hot object, growth of certain biological systems,and a host of other related first order phenomena. Despite the model’s simplicity, this kmd of problem occurs with surprising frequency in nature. If it is given a specified rate constant, the program yields amount versustime. If the input is normalized data, the program can be used to estimate a first order rate constant. As an illustration of its use (fig. 4), a graph is shown for thb ?utput of three different choices of the rate constant. The results are linear if plotted on semi-log paper. Some of the more interesting related examples might include the following: steady input (instsad of a step functianL periodic input. multivle first order inputs with different time constants, multiple first ord& leaks, a combinoton of multiple inputs and multiple I&s, a functionally dependent rate constant (espe-
.
.
.
The d&rential equations which apply for direct analysis, and their transforma!ion into incremental fom, and camputer language ?.re:
t.t=(-b-
1.1 + k-2*2,,)
= -k2,A,,, d4.r
=(kaAl,r
-3
dr + k,,A,,, - knAz,r)
dt . d’
=k2,A,,rdr-k12*2,,dl, Dl =-Kl
*A(I,I)“D4+#2’A(2,1)*D9,
D2=
Kl”A(I,I)‘“D4-K2*A[?,l)*D9,
Ai.r+t
=At.r+~t.r,
A2,,+,
=4.r
+ d4.r
3
A(I.I+I)=A(I,I)+Dl, A(Z,I+l)
= A(2.1) + D2.
4.2.2. Treatment of the variables To identify the amount in each compartment at a specified time, double subscript notation is used in this program. This permits generalizalimt of the made1 and the introduction of matrix methods, if desired. (In most instances. matrix methods and hiehcr mathec&w they might wve as a simplification). Initnlization is incorporated into the direct analysis calcolations. Values of the variables are computed and stored io place of prior values calculated from the direct analysis equation. The incrementation is accomplished by calculating and adding small changes in the direct analysis variables as they pass through the loops of the program. 4.3. Discussion The results are graphed as an cmnt in each compartment versus time, as shown :n fig. 5. If the rate cotwants are known. the amow.ts and/or concentrations in each compartment can be calculated. If the concentmtioos or amounts are known, an iterative procedure may be used to detennine the rate COD staots with the same program.
Several mterestingvsriations may be considered, ??afollows: I) he values and ratios of the rate constants may be varied sod the effect upon the cooteots of the compartments noted; 2) the rate eoo~tant~ may be replaced by functional relations; 3) these functional relations are permitted to interact: 4) the step function input may be replaced by another kind of input, such as continuous monotonically iowasing or decreasing (linear or non-linear) input; 5) leak paths may be added from the compartments to the exterior; 6) based on this program, a model may be developed which includes the following: a) passive diffusion between two compartments; b) active pornping betwan two compartments: c) leaks; d) “facilitated” diffusion; e) a coupling pump which rotates concenttation values and/or leak rates. In summary, this two-zomp~~tment method may represent such biological situations as the exchange of calcium between plasma and skeleton, the exchange of oxygen between cells and plasma. or the movement of a tracer or other substance between two body poo!s (as nitrogen between blood and muscle). Often, the two compartments are of unequal size, and the trans fer coeffcients are of unequal mayitode.
5. Compartment CQMP4)
analysis, part 3 (program
name:
5.1. Statement of the problem The distribution of body water in a human subject can be represented by a four-compartment model (Moore, 1962). It is desired to explore the model for various choices of mte constants. The diagram in fig. 6 shows the model chosen. 5.2. Fomr for direct analysis 7he physical laws in differential
form are:
The usual techniqueis followed in ewodhg the physicsl 1s~. The indcpadent whblc jt time, tad
tmountseach
the depsndcntvariabIes sre the In cmtprttnent. The model pammaters are the rateconstants,the volttmes, Qe twmber and
compttme11t
kind ofintereompttmentrl canttectionr, and the site ttttdnatureof the introductionof thedose of the tmtmialto be awed. The output it a table of the
103
FlNlTE ELEMLNT METHODS :N CELL DYNAMICS m the syrir
which may be fotknved by a finite
etenent analysis program. If one injects tritiated water into a subjjt and ML an iterative profedure to estimate r*te constants and compartment volumes, one&aim the valuesgiven in table 5. The data was obtained from normal human sub jcctr uoingradioarsays of samples which wem drawn serially from blood, urine, and pleural water. It can be shown that this modelcan assist differential discnwis ofdkeaw, guide and monitor therapy, evaluate health, and cixsSy diseasedstates.For example, in diabetes. urine outflow ir increased, so that et
is
matkcdly elevated;ln body dehydration, the volumr
t +A,aremall): during wrercise. there if atemporary shiftof w3tcr of eotnpattments L depressed(d
from eanpartmmt 1 to compnment 2;during heavy exertion, this ir accompmded by increawd breathing with cacommitlnt tisc ink3,; in uremia, kS1 ir dec~ea.wd; in edema 01 pqnancy, is added to be atialyzcd: in
amtp~cnf
cation. the volttmc of wmpwtnwd
another
walerint0.G
I increases.
6. JMueim 6.1. Gnuml~~r lxftwml md omto(isBI* traqwt proEctsn impottoot in bakgtcet systemr Tbeae proeosra, of tnetebolitee into ot out of e cell.
wdt
inthemovcwttt
If well b#ed in e hvpcrtonk or hypotonk scdw th. the wnwnlmtlon of coratituente in Ike interior of the all will be attcred bv diffir-ion and bv osnmsir
chatertiacmntbeceumnlh~*a.IM m&y
mthormlkd psrfomwd
lions For a&
analysi of dlffwton by solving parfii
yt of baundaty
is can.
diffarential
conditions,
qua-
m in-
diWtwl rdutkxt may be biwdaod concentration pmaeemaykeakeltedurfoaetionoftims (EkqerudI)lvidc. 1965). Tbeptob&mnmybeeppmedledbytiniteekmont lecbniquee (Conk. 1956). The division OFthe model into @tile danente is iilu~ted in fie 8. Our eartic-
arc connected in skea, aa iii In a chain.) The bowd&s between individual compattmmtr are ewned to be membranes pemdttmg sohrmt, but
notdutss, topass. (Animpermsblc bwnday nuy btplwd in Be model toihtntc the effwtoPk4kwing ncitticr alutenorsolvent topas.) ThS computmentc are identiiiid wtztNy by UK lubrript i, which ‘iantes from unity to i,,,, the aptetfud number dormpert-
mento. ‘MStypeof analysis
total
ml&
be useful for studying
prersure pulses traveling down an art&y.
Far tneny
bioh&ioel system8p~eteiq dunping, I ir pwibte toextend the analyeie by adding auitabledtunpimg factorr. In tlying to solve tbe diffusion prablm for timedcptndent Idutione, clo&aUy, one mult e&e the difffis;on equation !br the geometry md boundary conditions of intenet. When ooo trier to do thiefor multiplcccrmpartni6ntsmd multiple components, the rolutian bccomcc compltwed and difficult.
R.M0nltE,
Log
N.DAYIDS
Finite elemenl analysis simpliiics the conceptualization of lbe problem and the calcolatioo of answers, and one is not limited to the simple gccmetries. The solute diffuses across the semi-permcahle inarhces between the commirtments. Iniriall~~, it is all in compartment 1 6.2. The physicat Lws 6.2.1. Giffusion The physical law describing diffusion is Fick’s first law: the ner flux of scdute across an interface between two cell is proportmnal to the concentration gradient: Ji = -NC,+,
ci)Ar/Pr
,
where Ji is tbc amount of solute, in moles, trussing the right boandary of the itb cell in the positive direction per unit area in time Ahr(xc);IJ is the diffi;sion coefficient (cm2/sec),D is also the rate of flow scrouunit area for a unit difference ofconcentration per unit thickness:ci is the concentmtion of ith cell in mdes/cm3; Ax is the distance belween the midpoints of the two adjacent compartments. 6.2.2. Comeration cf mass The compartments are assumed to be finite in volume. Consequently, when wlute flows into, or out of, R compartment, the con.xntration in ‘bat cornpal Imeni is chnnged. The law of conservation of mass stae~ that the total mass in the system is constant. If mal~:rial flows out of one compartment, it flows into andhrr compartment. This law is employed to calculute new concentrations in the compartments xaording to the following mass balances: mass 10~ by flux out of the right boundary: massgain by flux in at left bouw dsry: net massgain: original mass in cell: new mass in Celi: net Ina*, gain:
equating
the two mass gains:
AFt AFt-1 -A(Ft @AX
Ft.. 1)
6.3. The analysis 6.3. I. Initialization Tbe values of the constants of tie medium, ,P), Ax. in,, and the value of Ar hie chosen. T;xe diflusion proce~ is then initiated hy assigning arbitrary starting values, c,, i = I, .,., im, to the concentration of the diffusing species in each compartment at time Since the left boundary of cell I ir assumed to be inpermeable, F, = 0. the following sequence of c:dcclations is performed:
I =0.
(I)
F, = -L&+,-
(2)
c;=q-(Fi-F,_,j/Ar,
(3)
replace ci by the new c; c&dated in step (2).
ci)Ar/oS
. i= l,._,im
If these calculations are carried out in the proper order, each quantity becomes defvled by a previous step. The last time these equations are used is when i= i,,,. The flux into the ([_ + 1)sl. cell is calcul&ted, but the concentmtion of that all is never provided an opportunity to chance. Hcoce it behaves as an infinite rcsewoir. After the amount Fi of the solute flowing across the ith interface has been cslculaled. the time is incremented by replacing the torrent time f by 116new value t’=l+Ar,
f=t’
and repeating the calculations (1)43). We continue in ibis manner for a specified number of time intervals. If a problem is sufficiently complicated, one may worry about instability in the solution and ti:.c appearance of undesired oscillations. This problem has been considered by Crank (1956). and he employs a diiensionless modulus which is a ratio between variables. If the value of this modulus is kept less than a certain value, the solution remains stable, and the values of the variables arc adjusted to satis) this condition.
CiALW (c; - q)AW
-A(F, -x=-t, =(c:-cij4AX This relation is used to solve for the new concentration. c;: c; = ci - (Fi - Ft_ +X’.
and R.I..BERGER
The above analysis is in a form which facililates writing a computer program (Dwids wd Berger, 1964). Except for notation. the ldtyskal laws, stated in (I), (2). and (3) above, are in 21form which can be used as statements in the progmm. almost in the fan in which they stand. The notation most be consistent
with the requirements of the particular programming system chosen.To carry out the sequenceof calculaLionsfor i = loop most be provided, which ;E, in tom, enclosed in an outer loop to repeat all the required stepsfor a given time interval. Thus, in the “BASIC” version given in the appendix, appropriate “FOR” and “NEXT- statements arc used to define theseloops. In addition, suit bb “READ” :statements are omvided to initialize the oroblem and “PRINT” statementsto display the output, as well as to print labels for column headings.(A “READ” ~fatemeat CZ.WP the computer to go to B “DATA” altemeot. and assignthe numbers found there to the valuesof the variables listed ic the “READ” statement. A “PRIYT” statement causesthe computer to print or type the current valuesof the variables which are listcd in the “PRINT” statement.) It is cwvenient also to provide in this, as well as each of the other pr& grams, suitable documentation in the form of a symbol table and a iabel for each physical law, so is to make the analysis self-identifying from the program listing. This may be done by including remarks(identified bv “REM”) which do not causeinstruciions to be executed but which remain in the prrr grm for explanation, amplification, and &&cation. Other boundary conditionsmay easily be inserted, as shown by examples in Je results. The output of this program is in rhe form of an array, in which each row givesthe concentration in each of the ten compartments at a particul~~rtime. The concentrations at later times are printed in SWcessiverows. In addition, the vahtesare nonnaltied to be numbers between 0 and IGQCIfor conv&ence in printing.
I, _.., i,,,a
6.5. Redts The output data may be plotted in any of three ways: I) concentration versustime for each compartment; 2) concentration (at a specitii time) verms compartment; 3) compartment m which concentration reachesa specitic value wrsus time for each conclmtration, or percentage of final concentration. The tirst of the above +ernate ways was ued for the curvesshown below. l&y might be obtained experimentally by withdrawing samplesfrom each compartment according to a selected schedule, and assayed
for concentration of soluteiin order ihat the sampling not disturb the system, one must assumethat the volume of solvent and amount of solute in the sample withdrawn is very much less than the original volume of solvent and rolute in a compartment. (It is possible, however, to take into account a change tn dimensions of compartments.) Fig. 8 show the concentration prclilesversus time for canpartments for the caseof an impermeable barrier on the left, and initial concentmtions as shown. Note the interesting result, most pronounced for compartment 2, lhat the valoe of the concemration reachesa maximom in 2 to IO. and ther! declines. A run was also made (not shown) -tith ths impermeable barrier on the right and initial conditions reversed.As expected, the diffusion takes place in the reversedirection, and the graph ;f the numerical values is similar to the graph of the previous case. Fig. 9 shows the concentration orofiies which result ifdiffuion is aww from cell ID m-both directions. Here we have an _ infinite reservoir next to the ceil 0:‘inifial high concentration.
I to8,
The curvesfor fig. 9 represent the case where the impermeable barrier has been placed between cells 5 and 6. I” thecomputer program, thiscondition is accomplishedby inserting the ~tatcment F(S) = 0 (no flux between thew zclls) immediatciy after the &tcmcnl drtining f$. This &traics :hd easewith whi:h chungesin the haundnry conditions cc” be mad: wins a finite cl ment analysis. A similar analysis might be applicable to biomedical problems. The rate of diffu don of a drug through a tiswc may be described by concentration profiles sinule~ to those &“m above. %e area A wooid be assignedthe value of the equivclent area of a biologic cell hr would be bssigcd the equivalet$ thickness of a cell, and an appropriate value would be assignedto the diffusion coeffiiicnt, D. AL an example. fur crythrocytes in a rouleau, we have A = 50 square micra Ar = 2 micro P w,er = 0.236 0stn01es/cm%ec Concentratioo-time raves would then indicate when the concentration of solute in a time leached a selected threshold value. A bone would constitute a diffusion harrier. Similar cutvcsmight also be useful io irlstmmentatior,. 8s in tbc desirn of an artiticial
Wave propagation occursin physiological processes. So& processesmay be i\lustrated by a pressurepulse trnveling dam a” artxy, an cbmpt shock transmitted to the skeleton, an acoustic signal in the auditory mcctus. and electrical phenomena associatedwith ncrvc conduction. We sliall deszcribc here c finite clnncnt analysis for such phenomena (Mehtn and Davids, 1966).
7.1. Finirr eiemenl analysis A one.dimc”siotql system is assumedto be divided into a &es of equal compartments. The coolpartmtrtts arc arrigocd numbers, from I throw& IL The
width of each compartment is Ax The compartment treated by the calculations is f’*signated es the ifh compartment. The volume of the ~th c”mpartmc”t is ui. There 1%a pressureexerted on the left side “1’the ith compartment, and it is denoted by P,, and ttere is e” opposmgpressureon the rigbc side of tbe ith compartment, denoted by P!,,. lllc cell cc” move cs a result of unbdanccd forces, and it may acquire a velocity. For the ith cell, this i? called oi, For stabili’y, Ax is made equal to cAf, with c = dB/,x Each element interacts with its neighbor, according to the followng ,;ysternof re!ations: I. Volomc change of cell: &a :=(“,+t - vi)A.z’Ax 2. Resulting ncwvokooe: “i+, = ui+t + &!A 3. Pressure-volumelaw: P,+t i: B&+, 4. Coundary condition: Pi&t = 0 5. Impulssmomentum conservation law for velocity change: & = (-Pi+, +Pi)At/Ax o. Resulting new velocity: oi = U$+ nu Here B is the hulk modulus, “I presmrc-voleme cociiicient,p is the pressure,u is the valume of the medium, and p is the density of the substance.
Little else is needed to develop a computer prcgram bared on the aboveanalyst. Tie is taken to be ttte indcpendenr variable. and the pressureand wbxity in each compartment are depndrnt variables. Loops are provided to scanthe above laws thro‘;gh lhe cells of the medium and to repeat the wan for successive tex intervals in the samemannet cs in the diffusion program. Documentation and a list of notations is provided in the listing in “BASIC” which is avzdlable from the authors on request. Table 6 shows a sample result of the wwe propagation piogrem for a unit step inpu,: on the left edge of the medium. Note that the boundary condition of zero pressureet the right end gcneratcsautonnticaliy a reflected wave from right to left. This wave then zig-zagsback and forth acrossthe medium as theory requires.
8. Heal conduction Heat conduction plays a fondcmental role in all control processesof living systems,since it gwt&
HI61 99
0
99
0.
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 0 0
.
99
99 99
99
0
0
99
99
99‘
0
0
99
99
99
99‘
Q
0
99
99
99
99
99
99
-0.
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
YY
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
9')
9!,
99
99
99
99
99
93
99
YY
99
99
99
99
99
99
99
99
99
99
9'3
99
99
99
99
99
99
99
99
99
99
99
99
9_s_,--6
Y?
_y? ' 0
99. _*0
99, .*0
99
99 99/--O
'-9
--a.
0
0
0
0
0
0
-0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
cl
0
99 ‘-..q 99
99.-.._o
99
99
9s‘~.,a
99
99
99
55.
0
Y?
99
99
99
33.
99
99
99
99
99
-;9.
99
99
99
99
99
99
.‘99.
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
'99
99
99
YY
99
99
99
'99
99
99
99
99
99
94
39
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
__O0
99
99
99
Q_=.
99
99
a..--'0
99
0,'
9_9_--
-0
99
0
-0,’ 0
,I
0
0
0
0
0
1
influences the operating tempere~re. A tinite element rnelysis permits addin this process to others which may be going on in I system at the same lime. Heat conduction is mathematically similar to diffusion, and the general differential equations describing each are similar (Crank, 1956; Car&w and Jaeger. 1959).
d
“,.“*.U3.U4.“..25*
-__
The model is assumed to be a semi-mfinite homegeneous slab with prescribed temperatures at the boundvies The bottom of SIIP!I i slab is to be immersec, in a lieuid nitrown bath while the top of the rleb is in rwm air. The temper: two profdes in the interior of the slab ver.sw he are d&led.
.
_
6.2. Finite element analysis The details of this analys:s are jive” by Berger and Davi&(1965), and Davidsand Berg01 (1964). Example A graph of output data (temperature versus lbne) may be found in fig. 10. Inhially the temperetcre of the bottom IO compartmenl~ wes taken to he that of the liquid nitrogen bath (- l96’C). the temperature of the top 5 compartments to be ambient temperature (+2W). Time is ndwnced in SO millisecond increments. This graph shows that the tern peratwe in comprttrent I and compartment 16, is taken to be at r~~rn temperature nnd liquid nitrogen tetnpereture, respectively. The intermediate campartments ultimately reach P linear strady state, varying betwen the boundary temperatures. The most interesting coml;artm@~ts a-e numbers 5 end 6, which initially support a tempemture difference of 221% This c?mot long cohtir.ue, and the curves show bow the magnitude of this gradient decays. This analysis is of interest in heat transfer pmb Iems, where temperature profiles ece Aeeded end freea ing rates are to be calculated to produce these. The free&g of blood platelets and the inf!ux of heag @xv.@ the wells of e storage vessel or shipping con-
the
tainer cooled by liquid nitrogen arc biological applications. with different edge conditiqls, the tempe!ature protIles might describe the dissipatiotl of heat from a working * tale. or, in reectar engineering, the loss of heat from a nuclear reactor fuel element. The campartments can be made nvxe or le~ numerow than 15. lltey IP ay be assigned varied value(for thermal
conductivity, specific heat, heat capacity, or Ihickners, 01 they may start with a different initial temperature distribution. Carslaw end Jaeger (1956) have t:onsid; wed the problem of stabili& of solution and accurec:~ of numerical values. They define: a mod&s, the s alue of which mu%t be kept smell to insure stability.
9. Summary A method is described)?elled finite element analysis. which facilitates obtaining umerical solutiwts for mathematical models. It is slu? able for eflicient computer implementation, end it is based directly on the relevant phytical law(s), and its results BII: interpretable in terms oft& variablei of the problem. The method in applied to the following problems: :) a weight felling in a gravitational field,
5) diffusion of a substance, 6) propagation 0r a WBVC, 7; beat conduction The method does not gw solutions in closed ,‘wm, and experience is necessaryin designingthe loops tu place the stepsin the proper sequence.If the problem involvescyclic behavior. a test for instabihty and choice of the correct increment for the indcpesdent variable must be mzde.
Acknawledgment We wish to thank Robert Friedcnberg fw hrs maay editorial suggestmnsduring the preparatnn of :he manuwipt. This investigation was support,:d m part by Public Health Service ResearchGrants ros. HE-10320(HEM) :lnd IR-OIHEI 1289-02. from the Nat,onal Heart Inst,tute.
There approaches are in contrast to tbc classical methods of (a) mathematical analysis and(b) convcntional computer programming. Commonly, the laws governing a system are analyzed mlthematically by dcrtving differential equations and by solving them in closed form, wherever possible. These differential equations, or particular solutions of them, are ttad at the beginnir.g 0; suctt a mathematical analysts. In contrast, finite
R
MOORE. N.DAYIDS
and RLBERCER
studied. This may be accomplished by using 8 running index t” designate which parts oft problem arc being considered. The physical law, in incremcnta’l form, is tberl staed. Vsrious examples arc pen further in the text. Tbc initialization of the variahlrs is necessary t” rtiot tbc cycle ofcom~utations, and a set of bandary conditions is needed which arc similar to those in clas:ical mathematical analysis. In the example of a freely falling body, to be discussed. if a weir&t has ari initial velocitv of zero. and if it falls from a heaht of one hundred feet, the”: u= O,andy = Ior,, for r=O. An important part of the calculations ia fimte element analysis involves incrementing the variables appearing in the physical laws and summing over these increments by looping back in the program. The magnitudes of the increments arc cb6s:n to satisfy criteria of soccitied accuracy and “erhaps #criteria of stability. In a-given problcm,wrt& ratios, referred to a~ “moduli”, arc useful in f&ing the ratio “f the incremcots, such as that between distance end time. The next step in the constmctit~n of a finite ~clcment analysis program is to establish the loops by which the summation over the vaoables is performed. A “loop” consists of ““c or m”rc
xe
.
F&z. i. scieme for problem soI”¶iO” by ahnate
methodi.
In conventional approaches to the solution oim:tthcmatical forms of biophysical problems, the invcstigator may start from a previously derived mathcmatic.d solution in closed form. expressed a an mtegr;d “I as a series. Hc may ultimately apply a computer to helo him. A “finite difCeerence” method. or other nunwical method, is frequently used to arnve at a solution. w!wn such a closed-form solution is unavailable. Finite clement analysis, by cootmst, operates on tbc laws themselves, not on derived relations. The following steps arc necessary in carrying out the finite element attalysis of a pm5lem: i) statcmcnt of the probLcm, 2) subdivision of the physical structure into finite elcmcna, 3) initialization of the variables, 4) calculation in sequence (by “loops”), a.. exprcssirl of the physical laws in incremental form. b. the summa!ion of the wriablx The statemen: of the problem mus? be phrased in ss specific terms as possible. Thus, the transport process thraugh cellular membranes may bc resolved into a :argc number of specific components. It is useful to desigr.ate the specific components which arc being
3. Averagingof initial and final velocity over the time interval: ” = ;(ui + q). 4. (a) Incrementing the height: Ayv=uAt. (b) Accumulation of height mcrementr, yf =yi + Ay. 5. We then mcrement tho time, 1, by a step. AI, and repeat the above cycle of computations. In this simple case, the force is assumedto be constant. If the force is not constant, as in simple harmonic motion, then correction terms have to be used to calculate the final velocity. The discussionof this problem occurs ekewhere (Balko and Berger, 1969). A suitable magnitude for the time increment can be establishedby trial. The problem is mn with EUCcessivelysmaller time increments, until the graphs of successivesolutions superimpose.Once the solution is found to be invariant with respect to choice of time increment, the value for the time increment may bs mcreasedto conservecomputer time, but still to obtain a solution within the desired bounds of accuracy. Thts pant will be dealt with in depth in later papers. The solution to a simple problem is: A” =gAf = 32.2 X 0.1 = -3.22 ft/sec , uf = u, + Au = 0 - 3.22 = -3.22 It/xc u=t(ut+u,)=:(O-3.22)=-I.61 2. Simple application
Ay==uAr=-1.61
X0.1 =-0.161
, ftlsec,
ft,
yf =yi + Ay= IO0 - 0.161 = 99.R36 ft. This section will describe the application of the finite element approach to a simple problem, the mu tion of a freely-falling body. Later (Moore, Gwids and Berger, 1969; Balk0 and 13erger,1969), more complex problems, of interest to biology, will be dedt with, since the greatest advantageof the method is its apphability to highly complex systems. 2.1. Andysis of mor;oon in c gmvimtioml fidd :. Initialization ofvariables: At the beginnmg softhe problem, the time, I = 0, and the velocti~i. ” = 0, and the = IOQ feet. 2. Incrementing the variables: (a) the velocity iwement, aV =gAt, (b) the accumulation of velocity is done by addmg the increment, Au, to the initial value ~oCvelocity. u,, to get the final value of velocity. q: uf=ui+Av.
height,y
rf=ti+Ar=O+O.l=O.lsec. Ii the redder cares to perform this set of calculations for twenty-five times, he will obtain a table of results similar to the output of the computer program. 2.2. Computerpwgmnt The jldzated cumulation of the variables may be performed by machine. This is done by including these operations in a loop. A loop is a set of instrw Lionswhich are executed one or more times. The cumber of times the calculations in the loop are repeated may be controlled in severalways. This section will discussthe kinds of “BASIC” statements that may be used to control a loop: 1) “FOILNEXT” statements 2) “LF...THEN” state!ne”ts 3) Independent mop counter
R.MOORE.N.DA”IDS
98
A “FOR...NEXT”’
statement is a pair of instmc-
tions of fhc form: sI=IT03
STEP I
I i!i tbc indcpcndent variable; the “FOR” statement MS its value to I, and the subsequentinstructions are perFormed down to the “NEXT I” statement. The program goesback to the “FOR” statement. increases the xdue of the I by i from the instructions to the ‘NEXT” statement, goesback to the “FOR” slatemer.t, and increasestbc value ofI by from 2 to 3, at which point I equals the second nunbcr in the “FOR” statement. The instructions arc performed for the last time. curdthe prcaam continues prist the “NEXT I” statement, to the remainder of the pmgratn. In the pr5grun cntitlcd “GRVTY I”, “FOR...NFXT” statanents ue used in lines 520 and 6 IO. Fig. :#showsa conpanson between the analytic and the finite slemen: solutions.
I tr,2,exewtes
1
and R.L.BBRGER
condition to tc met. If the condition is met, then the prol~zm ;s directed to a line, identified by a number, to cw tinue the calculation. If the con&i!on is not met, then the next ins:ruction in !equencc :s executed. For c:xnmplc,“IF T < 2.6 THEN .jZo” will cumparc the a rrcnt value of T with 2.6. If T is less than 2.6. then the comptttcr goesback to Ii le 520 and ~‘xewtcs that tolstmction.1n”GRVTY 2”,“IF...THEN” stcbcmew! are used in lines 520 and 6’10. In “GRVTY 3”. an independent loop counter is employed, rather than looping ba :cd an the valuesof any OFthe val-iables,as in lines 520 and 590. It should he noted that the symbols Dl,DZ,... are intendsd to stand for At nnd Ax. becauseof the unavailability oi the Greek “A” symbol and the limitations in naming I variable in the “BASIC” langual;c. However, D I. D2,.. should not be confused with differentials. The sectitmsof the program arc mtmb’eredand labeled in which major portions of the work ore done. The subdivisions of the computations section, by groups of instructions, are shown in the table below.
F%yri:allaw: Initkilting tll~ “1riable-. tllCt~W”ting tie “ti*bt:l Lwp controt
The solution was checked by calculating the six _ displscemcntsfrom tie cxact formula y = 100 -OS&. f = 0.1,0.2, .. .. 2.5. The results agreed to three decimal placer. An “IF...THBN...” statement is one which states a
530.s70
53O,S?O
550.580
310-370
3to-370
310.-3,0
540,560 520.610
520.540.580 520,6,0
540.S60.590 5ZO.610
III other words, the three gravity programs illustrate due: methods of looping. In I, the looping is controlled by the instruction, “520 FORT = 0.1 TO 2.5 STEP Dl”. This setsthe variable T to an initial value of0.l. executes the instructions down to “610 NEXT T”, and returns to line 520. The value of T is increasedby the amount DI (“331) LET 01 = 0.1’3, ar-i the instructions are again executed to line 610. Thi:. is done until T= 2.5, for which the instructions arc xccuted fur the final time. In II, the control of the loop is different. and it is done by”610 IFTC 2.6THEN 52O”.Tbe initial ~~al~eofTI;(T=T+DI=O+O.l)=O.i.Theinsrructior s are executed down to line 610. This compares
the value of T wth 2.6. If the value of T is less than 2.6, the program returns to line 520. If T is equal lo 3.6, then the program stops. In 111, the looping is controlled by line “530 FOR I = 0 TO 24”. The variables HPCsubscripted, and the value: of the subscript, I. increases by 1 for rach time incren~t, until the value of the subscript equals 24 Then the instructions are executed for the linal time. I is, in effect, a loop counter which keeps tlack of the number of times the program goes thm the set of calculations in lmes 540-6GO. There are additional kinds ofcontrol of lowine possible. One may have several nested loop:;. l&g combinations of the above three kinds of control. and there may be severa: nested loop counters. There may be comparison statements for variables othsr than the independent variable. T&se may be used so as to causr the program to generate families of solutio~a. for many ranges of the variables, and even to search for, and to print out, only those few solutions which are of interest for some particular reason. The availability of several means of conrrollmg loops frees the investigator from limitations inherent in either the notation employed or the mbrcnpts allowed by the language.
3 Compartmental
where 01 = change of amount of substance in the corn. partmen,, IL= rateconstant for the process, A(1) = amount of substance in the comoartmrnt at the Itb time. D2 = time mcrement. After DI has been caiculated. the increment ofA is computed, and this increment is added to the previous value A. This new value computed for A is ~tilbed m the next increment reque&This is accomplished by the instruction “LET A(l+I)=AlIPDl”. ‘Tlme is then incremented, using the in&c& “LET T(l+ I) =T(I)+DZ”. This also performs the operation of loop ing back for successive increments. The program then calculates A(0 versw r(r) for three different values of the first order rate constant. k, as shown in fig. 4.
of
analysis, part 1
Compartmental analysis may be used to describe the movement or change of substances in luantitattve fashion (Hart. 1963). For example. finite clement analysis enables the reduction of data from radioactive tracer studies to be performed in a very simple WY.
In a chemical :eaction in whxh A gces to B, v&ten A&B, the rate equation governing this process, called first order kinetics, is given by &I = -kAAi
:.te
where bA is the change in amount, A, in cornpartment in the time At; k is the first order rate constaM inset-t. In “BASIC” notation this is writ’en as: DI = -l?A(I)*DZ,
The variables are handled in the following way: time is considered the independent variable, and the amount of substance present in a cotr.partment is the dependent variable. The amount initially present, and the rate constant, are program parameters which may be changed to fit various experimental conditions. 3.3. Dmmenrarion The documentation perts:
consists of the following three
cially one which is B function of :he input rate), and coup!ing between lass rate and input rate by mwlr of wriox functions (with and witbout time delays). Each casxcan be explored by changing just ora: to three lines of the program.
4. Compartmental analysis. part :! (prugmm name: COMP2)
:;tatemevr ofrhe
4.1. problem The problem is to determine the kinetic bebaviaur of a ‘we-compartment system which cxchaogescontents by first order Proresser(see fig. 5).
3.4. Discussion This case is simple but common. It describesrsdi* isotopic decay, first order chemical reactions, operation of a shock absorber, cooling of a hot object, growth of certain biological systems,and a host of other related first order phenomena. Despite the model’s simplicity, this kmd of problem occurs with surprising frequency in nature. If it is given a specified rate constant, the program yields amount versustime. If the input is normalized data, the program can be used to estimate a first order rate constant. As an illustration of its use (fig. 4), a graph is shown for thb ?utput of three different choices of the rate constant. The results are linear if plotted on semi-log paper. Some of the more interesting related examples might include the following: steady input (instsad of a step functianL periodic input. multivle first order inputs with different time constants, multiple first ord& leaks, a combinoton of multiple inputs and multiple I&s, a functionally dependent rate constant (espe-
.
.
.
The d&rential equations which apply for direct analysis, and their transforma!ion into incremental fom, and camputer language ?.re:
t.t=(-b-
1.1 + k-2*2,,)
= -k2,A,,, d4.r
=(kaAl,r
-3
dr + k,,A,,, - knAz,r)
dt . d’
=k2,A,,rdr-k12*2,,dl, Dl =-Kl
*A(I,I)“D4+#2’A(2,1)*D9,
D2=
Kl”A(I,I)‘“D4-K2*A[?,l)*D9,
Ai.r+t
=At.r+~t.r,
A2,,+,
=4.r
+ d4.r
3
A(I.I+I)=A(I,I)+Dl, A(Z,I+l)
= A(2.1) + D2.
4.2.2. Treatment of the variables To identify the amount in each compartment at a specified time, double subscript notation is used in this program. This permits generalizalimt of the made1 and the introduction of matrix methods, if desired. (In most instances. matrix methods and hiehcr mathec&w they might wve as a simplification). Initnlization is incorporated into the direct analysis calcolations. Values of the variables are computed and stored io place of prior values calculated from the direct analysis equation. The incrementation is accomplished by calculating and adding small changes in the direct analysis variables as they pass through the loops of the program. 4.3. Discussion The results are graphed as an cmnt in each compartment versus time, as shown :n fig. 5. If the rate cotwants are known. the amow.ts and/or concentrations in each compartment can be calculated. If the concentmtioos or amounts are known, an iterative procedure may be used to detennine the rate COD staots with the same program.
Several mterestingvsriations may be considered, ??afollows: I) he values and ratios of the rate constants may be varied sod the effect upon the cooteots of the compartments noted; 2) the rate eoo~tant~ may be replaced by functional relations; 3) these functional relations are permitted to interact: 4) the step function input may be replaced by another kind of input, such as continuous monotonically iowasing or decreasing (linear or non-linear) input; 5) leak paths may be added from the compartments to the exterior; 6) based on this program, a model may be developed which includes the following: a) passive diffusion between two compartments; b) active pornping betwan two compartments: c) leaks; d) “facilitated” diffusion; e) a coupling pump which rotates concenttation values and/or leak rates. In summary, this two-zomp~~tment method may represent such biological situations as the exchange of calcium between plasma and skeleton, the exchange of oxygen between cells and plasma. or the movement of a tracer or other substance between two body poo!s (as nitrogen between blood and muscle). Often, the two compartments are of unequal size, and the trans fer coeffcients are of unequal mayitode.
5. Compartment CQMP4)
analysis, part 3 (program
name:
5.1. Statement of the problem The distribution of body water in a human subject can be represented by a four-compartment model (Moore, 1962). It is desired to explore the model for various choices of mte constants. The diagram in fig. 6 shows the model chosen. 5.2. Fomr for direct analysis 7he physical laws in differential
form are:
The usual techniqueis followed in ewodhg the physicsl 1s~. The indcpadent whblc jt time, tad
tmountseach
the depsndcntvariabIes sre the In cmtprttnent. The model pammaters are the rateconstants,the volttmes, Qe twmber and
compttme11t
kind ofintereompttmentrl canttectionr, and the site ttttdnatureof the introductionof thedose of the tmtmialto be awed. The output it a table of the
103
FlNlTE ELEMLNT METHODS :N CELL DYNAMICS m the syrir
which may be fotknved by a finite
etenent analysis program. If one injects tritiated water into a subjjt and ML an iterative profedure to estimate r*te constants and compartment volumes, one&aim the valuesgiven in table 5. The data was obtained from normal human sub jcctr uoingradioarsays of samples which wem drawn serially from blood, urine, and pleural water. It can be shown that this modelcan assist differential discnwis ofdkeaw, guide and monitor therapy, evaluate health, and cixsSy diseasedstates.For example, in diabetes. urine outflow ir increased, so that et
is
matkcdly elevated;ln body dehydration, the volumr
t +A,aremall): during wrercise. there if atemporary shiftof w3tcr of eotnpattments L depressed(d
from eanpartmmt 1 to compnment 2;during heavy exertion, this ir accompmded by increawd breathing with cacommitlnt tisc ink3,; in uremia, kS1 ir dec~ea.wd; in edema 01 pqnancy, is added to be atialyzcd: in
amtp~cnf
cation. the volttmc of wmpwtnwd
another
walerint0.G
I increases.
6. JMueim 6.1. Gnuml~~r lxftwml md omto(isBI* traqwt proEctsn impottoot in bakgtcet systemr Tbeae proeosra, of tnetebolitee into ot out of e cell.
wdt
inthemovcwttt
If well b#ed in e hvpcrtonk or hypotonk scdw th. the wnwnlmtlon of coratituente in Ike interior of the all will be attcred bv diffir-ion and bv osnmsir
chatertiacmntbeceumnlh~*a.IM m&y
mthormlkd psrfomwd
lions For a&
analysi of dlffwton by solving parfii
yt of baundaty
is can.
diffarential
conditions,
qua-
m in-
diWtwl rdutkxt may be biwdaod concentration pmaeemaykeakeltedurfoaetionoftims (EkqerudI)lvidc. 1965). Tbeptob&mnmybeeppmedledbytiniteekmont lecbniquee (Conk. 1956). The division OFthe model into @tile danente is iilu~ted in fie 8. Our eartic-
arc connected in skea, aa iii In a chain.) The bowd&s between individual compattmmtr are ewned to be membranes pemdttmg sohrmt, but
notdutss, topass. (Animpermsblc bwnday nuy btplwd in Be model toihtntc the effwtoPk4kwing ncitticr alutenorsolvent topas.) ThS computmentc are identiiiid wtztNy by UK lubrript i, which ‘iantes from unity to i,,,, the aptetfud number dormpert-
mento. ‘MStypeof analysis
total
ml&
be useful for studying
prersure pulses traveling down an art&y.
Far tneny
bioh&ioel system8p~eteiq dunping, I ir pwibte toextend the analyeie by adding auitabledtunpimg factorr. In tlying to solve tbe diffusion prablm for timedcptndent Idutione, clo&aUy, one mult e&e the difffis;on equation !br the geometry md boundary conditions of intenet. When ooo trier to do thiefor multiplcccrmpartni6ntsmd multiple components, the rolutian bccomcc compltwed and difficult.
R.M0nltE,
Log
N.DAYIDS
Finite elemenl analysis simpliiics the conceptualization of lbe problem and the calcolatioo of answers, and one is not limited to the simple gccmetries. The solute diffuses across the semi-permcahle inarhces between the commirtments. Iniriall~~, it is all in compartment 1 6.2. The physicat Lws 6.2.1. Giffusion The physical law describing diffusion is Fick’s first law: the ner flux of scdute across an interface between two cell is proportmnal to the concentration gradient: Ji = -NC,+,
ci)Ar/Pr
,
where Ji is tbc amount of solute, in moles, trussing the right boandary of the itb cell in the positive direction per unit area in time Ahr(xc);IJ is the diffi;sion coefficient (cm2/sec),D is also the rate of flow scrouunit area for a unit difference ofconcentration per unit thickness:ci is the concentmtion of ith cell in mdes/cm3; Ax is the distance belween the midpoints of the two adjacent compartments. 6.2.2. Comeration cf mass The compartments are assumed to be finite in volume. Consequently, when wlute flows into, or out of, R compartment, the con.xntration in ‘bat cornpal Imeni is chnnged. The law of conservation of mass stae~ that the total mass in the system is constant. If mal~:rial flows out of one compartment, it flows into andhrr compartment. This law is employed to calculute new concentrations in the compartments xaording to the following mass balances: mass 10~ by flux out of the right boundary: massgain by flux in at left bouw dsry: net massgain: original mass in cell: new mass in Celi: net Ina*, gain:
equating
the two mass gains:
AFt AFt-1 -A(Ft @AX
Ft.. 1)
6.3. The analysis 6.3. I. Initialization Tbe values of the constants of tie medium, ,P), Ax. in,, and the value of Ar hie chosen. T;xe diflusion proce~ is then initiated hy assigning arbitrary starting values, c,, i = I, .,., im, to the concentration of the diffusing species in each compartment at time Since the left boundary of cell I ir assumed to be inpermeable, F, = 0. the following sequence of c:dcclations is performed:
I =0.
(I)
F, = -L&+,-
(2)
c;=q-(Fi-F,_,j/Ar,
(3)
replace ci by the new c; c&dated in step (2).
ci)Ar/oS
. i= l,._,im
If these calculations are carried out in the proper order, each quantity becomes defvled by a previous step. The last time these equations are used is when i= i,,,. The flux into the ([_ + 1)sl. cell is calcul&ted, but the concentmtion of that all is never provided an opportunity to chance. Hcoce it behaves as an infinite rcsewoir. After the amount Fi of the solute flowing across the ith interface has been cslculaled. the time is incremented by replacing the torrent time f by 116new value t’=l+Ar,
f=t’
and repeating the calculations (1)43). We continue in ibis manner for a specified number of time intervals. If a problem is sufficiently complicated, one may worry about instability in the solution and ti:.c appearance of undesired oscillations. This problem has been considered by Crank (1956). and he employs a diiensionless modulus which is a ratio between variables. If the value of this modulus is kept less than a certain value, the solution remains stable, and the values of the variables arc adjusted to satis) this condition.
CiALW (c; - q)AW
-A(F, -x=-t, =(c:-cij4AX This relation is used to solve for the new concentration. c;: c; = ci - (Fi - Ft_ +X’.
and R.I..BERGER
The above analysis is in a form which facililates writing a computer program (Dwids wd Berger, 1964). Except for notation. the ldtyskal laws, stated in (I), (2). and (3) above, are in 21form which can be used as statements in the progmm. almost in the fan in which they stand. The notation most be consistent
with the requirements of the particular programming system chosen.To carry out the sequenceof calculaLionsfor i = loop most be provided, which ;E, in tom, enclosed in an outer loop to repeat all the required stepsfor a given time interval. Thus, in the “BASIC” version given in the appendix, appropriate “FOR” and “NEXT- statements arc used to define theseloops. In addition, suit bb “READ” :statements are omvided to initialize the oroblem and “PRINT” statementsto display the output, as well as to print labels for column headings.(A “READ” ~fatemeat CZ.WP the computer to go to B “DATA” altemeot. and assignthe numbers found there to the valuesof the variables listed ic the “READ” statement. A “PRIYT” statement causesthe computer to print or type the current valuesof the variables which are listcd in the “PRINT” statement.) It is cwvenient also to provide in this, as well as each of the other pr& grams, suitable documentation in the form of a symbol table and a iabel for each physical law, so is to make the analysis self-identifying from the program listing. This may be done by including remarks(identified bv “REM”) which do not causeinstruciions to be executed but which remain in the prrr grm for explanation, amplification, and &&cation. Other boundary conditionsmay easily be inserted, as shown by examples in Je results. The output of this program is in rhe form of an array, in which each row givesthe concentration in each of the ten compartments at a particul~~rtime. The concentrations at later times are printed in SWcessiverows. In addition, the vahtesare nonnaltied to be numbers between 0 and IGQCIfor conv&ence in printing.
I, _.., i,,,a
6.5. Redts The output data may be plotted in any of three ways: I) concentration versustime for each compartment; 2) concentration (at a specitii time) verms compartment; 3) compartment m which concentration reachesa specitic value wrsus time for each conclmtration, or percentage of final concentration. The tirst of the above +ernate ways was ued for the curvesshown below. l&y might be obtained experimentally by withdrawing samplesfrom each compartment according to a selected schedule, and assayed
for concentration of soluteiin order ihat the sampling not disturb the system, one must assumethat the volume of solvent and amount of solute in the sample withdrawn is very much less than the original volume of solvent and rolute in a compartment. (It is possible, however, to take into account a change tn dimensions of compartments.) Fig. 8 show the concentration prclilesversus time for canpartments for the caseof an impermeable barrier on the left, and initial concentmtions as shown. Note the interesting result, most pronounced for compartment 2, lhat the valoe of the concemration reachesa maximom in 2 to IO. and ther! declines. A run was also made (not shown) -tith ths impermeable barrier on the right and initial conditions reversed.As expected, the diffusion takes place in the reversedirection, and the graph ;f the numerical values is similar to the graph of the previous case. Fig. 9 shows the concentration orofiies which result ifdiffuion is aww from cell ID m-both directions. Here we have an _ infinite reservoir next to the ceil 0:‘inifial high concentration.
I to8,
The curvesfor fig. 9 represent the case where the impermeable barrier has been placed between cells 5 and 6. I” thecomputer program, thiscondition is accomplishedby inserting the ~tatcment F(S) = 0 (no flux between thew zclls) immediatciy after the &tcmcnl drtining f$. This &traics :hd easewith whi:h chungesin the haundnry conditions cc” be mad: wins a finite cl ment analysis. A similar analysis might be applicable to biomedical problems. The rate of diffu don of a drug through a tiswc may be described by concentration profiles sinule~ to those &“m above. %e area A wooid be assignedthe value of the equivclent area of a biologic cell hr would be bssigcd the equivalet$ thickness of a cell, and an appropriate value would be assignedto the diffusion coeffiiicnt, D. AL an example. fur crythrocytes in a rouleau, we have A = 50 square micra Ar = 2 micro P w,er = 0.236 0stn01es/cm%ec Concentratioo-time raves would then indicate when the concentration of solute in a time leached a selected threshold value. A bone would constitute a diffusion harrier. Similar cutvcsmight also be useful io irlstmmentatior,. 8s in tbc desirn of an artiticial
Wave propagation occursin physiological processes. So& processesmay be i\lustrated by a pressurepulse trnveling dam a” artxy, an cbmpt shock transmitted to the skeleton, an acoustic signal in the auditory mcctus. and electrical phenomena associatedwith ncrvc conduction. We sliall deszcribc here c finite clnncnt analysis for such phenomena (Mehtn and Davids, 1966).
7.1. Finirr eiemenl analysis A one.dimc”siotql system is assumedto be divided into a &es of equal compartments. The coolpartmtrtts arc arrigocd numbers, from I throw& IL The
width of each compartment is Ax The compartment treated by the calculations is f’*signated es the ifh compartment. The volume of the ~th c”mpartmc”t is ui. There 1%a pressureexerted on the left side “1’the ith compartment, and it is denoted by P,, and ttere is e” opposmgpressureon the rigbc side of tbe ith compartment, denoted by P!,,. lllc cell cc” move cs a result of unbdanccd forces, and it may acquire a velocity. For the ith cell, this i? called oi, For stabili’y, Ax is made equal to cAf, with c = dB/,x Each element interacts with its neighbor, according to the followng ,;ysternof re!ations: I. Volomc change of cell: &a :=(“,+t - vi)A.z’Ax 2. Resulting ncwvokooe: “i+, = ui+t + &!A 3. Pressure-volumelaw: P,+t i: B&+, 4. Coundary condition: Pi&t = 0 5. Impulssmomentum conservation law for velocity change: & = (-Pi+, +Pi)At/Ax o. Resulting new velocity: oi = U$+ nu Here B is the hulk modulus, “I presmrc-voleme cociiicient,p is the pressure,u is the valume of the medium, and p is the density of the substance.
Little else is needed to develop a computer prcgram bared on the aboveanalyst. Tie is taken to be ttte indcpendenr variable. and the pressureand wbxity in each compartment are depndrnt variables. Loops are provided to scanthe above laws thro‘;gh lhe cells of the medium and to repeat the wan for successive tex intervals in the samemannet cs in the diffusion program. Documentation and a list of notations is provided in the listing in “BASIC” which is avzdlable from the authors on request. Table 6 shows a sample result of the wwe propagation piogrem for a unit step inpu,: on the left edge of the medium. Note that the boundary condition of zero pressureet the right end gcneratcsautonnticaliy a reflected wave from right to left. This wave then zig-zagsback and forth acrossthe medium as theory requires.
8. Heal conduction Heat conduction plays a fondcmental role in all control processesof living systems,since it gwt&
HI61 99
0
99
0.
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 0 0
.
99
99 99
99
0
0
99
99
99‘
0
0
99
99
99
99‘
Q
0
99
99
99
99
99
99
-0.
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
YY
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
9')
9!,
99
99
99
99
99
93
99
YY
99
99
99
99
99
99
99
99
99
99
9'3
99
99
99
99
99
99
99
99
99
99
99
99
9_s_,--6
Y?
_y? ' 0
99. _*0
99, .*0
99
99 99/--O
'-9
--a.
0
0
0
0
0
0
-0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
cl
0
99 ‘-..q 99
99.-.._o
99
99
9s‘~.,a
99
99
99
55.
0
Y?
99
99
99
33.
99
99
99
99
99
-;9.
99
99
99
99
99
99
.‘99.
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
'99
99
99
YY
99
99
99
'99
99
99
99
99
99
94
39
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
__O0
99
99
99
Q_=.
99
99
a..--'0
99
0,'
9_9_--
-0
99
0
-0,’ 0
,I
0
0
0
0
0
1
influences the operating tempere~re. A tinite element rnelysis permits addin this process to others which may be going on in I system at the same lime. Heat conduction is mathematically similar to diffusion, and the general differential equations describing each are similar (Crank, 1956; Car&w and Jaeger. 1959).
d
“,.“*.U3.U4.“..25*
-__
The model is assumed to be a semi-mfinite homegeneous slab with prescribed temperatures at the boundvies The bottom of SIIP!I i slab is to be immersec, in a lieuid nitrown bath while the top of the rleb is in rwm air. The temper: two profdes in the interior of the slab ver.sw he are d&led.
.
_
6.2. Finite element analysis The details of this analys:s are jive” by Berger and Davi&(1965), and Davidsand Berg01 (1964). Example A graph of output data (temperature versus lbne) may be found in fig. 10. Inhially the temperetcre of the bottom IO compartmenl~ wes taken to he that of the liquid nitrogen bath (- l96’C). the temperature of the top 5 compartments to be ambient temperature (+2W). Time is ndwnced in SO millisecond increments. This graph shows that the tern peratwe in comprttrent I and compartment 16, is taken to be at r~~rn temperature nnd liquid nitrogen tetnpereture, respectively. The intermediate campartments ultimately reach P linear strady state, varying betwen the boundary temperatures. The most interesting coml;artm@~ts a-e numbers 5 end 6, which initially support a tempemture difference of 221% This c?mot long cohtir.ue, and the curves show bow the magnitude of this gradient decays. This analysis is of interest in heat transfer pmb Iems, where temperature profiles ece Aeeded end freea ing rates are to be calculated to produce these. The free&g of blood platelets and the inf!ux of heag @xv.@ the wells of e storage vessel or shipping con-
the
tainer cooled by liquid nitrogen arc biological applications. with different edge conditiqls, the tempe!ature protIles might describe the dissipatiotl of heat from a working * tale. or, in reectar engineering, the loss of heat from a nuclear reactor fuel element. The campartments can be made nvxe or le~ numerow than 15. lltey IP ay be assigned varied value(for thermal
conductivity, specific heat, heat capacity, or Ihickners, 01 they may start with a different initial temperature distribution. Carslaw end Jaeger (1956) have t:onsid; wed the problem of stabili& of solution and accurec:~ of numerical values. They define: a mod&s, the s alue of which mu%t be kept smell to insure stability.
9. Summary A method is described)?elled finite element analysis. which facilitates obtaining umerical solutiwts for mathematical models. It is slu? able for eflicient computer implementation, end it is based directly on the relevant phytical law(s), and its results BII: interpretable in terms oft& variablei of the problem. The method in applied to the following problems: :) a weight felling in a gravitational field,
5) diffusion of a substance, 6) propagation 0r a WBVC, 7; beat conduction The method does not gw solutions in closed ,‘wm, and experience is necessaryin designingthe loops tu place the stepsin the proper sequence.If the problem involvescyclic behavior. a test for instabihty and choice of the correct increment for the indcpesdent variable must be mzde.
Acknawledgment We wish to thank Robert Friedcnberg fw hrs maay editorial suggestmnsduring the preparatnn of :he manuwipt. This investigation was support,:d m part by Public Health Service ResearchGrants ros. HE-10320(HEM) :lnd IR-OIHEI 1289-02. from the Nat,onal Heart Inst,tute.