Lecture 5. - Part 2
Compartmental Problems
by Reinaldo Baretti Machνn
www.geocities.com/serienumerica4
www.geocities.com/serienumerica3
www.geocities.com/serienumerica2
www.geocities.com/serienumerica
For questions and suggestions write to:
Reference :
by William Simon ; http://www.amazon.com/Mathematical-Techniques-Biology-Medicine-William/dp/0486652475/ref=sr_1_1?ie=UTF8&s=books&qid=1223731765&sr=8-1
2. http://geocities.com/serienumerica4/Lecture1Biolologymedicine.doc
3. Mathematical Biology: I. An Introduction (Interdisciplinary Applied Mathematics) by James D. Murray
4. http://www.geocities.com/serienumerica4/Lecture3Biolologymedicine.doc
5. Schaum's Outline of Advanced Mathematics for Engineers and Scientists by Murray R Spiegel
6. Schaum's Outline of Laplace Transforms by Murray R Spiegel
Examples and problems are taken from:
Chapter 5 of Reference 1. The problem numbers refer to the textbook listing.
Sections 7 and 8.
The following example uses arbitrary parameters for I , F ,C not reflecting physiological values.
In a one compartment problem with out flow F ~ liters/minute the initial concentration is given as C(0) = 1 mole/liter with an input
I(t) = 10 ~ moles/min 0 ≤ t ≤ 10 minute
I(t) = 0 t >10 minute
Write the DE and solve it. In an interval ∆t the gain and loss in the one compartment are
∆Q = (∆t) ( I F C ) , divide by V and take the limit ∆t→ 0 ,
lim (∆Q/(V ∆t)) = dC/dt = I(t)/V (F/V) C(t) .
Using finite differences we get the solution
C(tn) = C(tn-1) + (∆t) { I(tn-1)/V (F/V) C(tn-1) }
Data
V=8. liters ,F= 1.25liters/minute ,Input= 10. moles/minute
Fig. C(t) when the input is a step function of 10 minutes duration
FORTRAN code
c from page 108 MATh Tech Biol & Med.
real input
data v,input,tf /8.,10,20./
f=10./v
tscale=V/F
dt=tscale/1000.
nstep=int(tf/dt)
print*,'nstep', nstep
kp=int(float(nstep)/60.)
kount=kp
c0=1.
Print*,'V,F,Input=', V,F,Input
print*,' '
print 100 ,0. ,c0 ,input
do 10 i=1,nstep
t=dt*float(i)
if(t.gt.10.)input=0.
c1=c0+dt*(input/V-(F/V)*c0)
if(i.eq.kount)then
print 100 ,t ,c1 ,input
kount=kount+kp
endif
c0=c1
10 continue
100 format(1x,'t,c(moles/liter),Input=',3(3x,e10.3))
stop
end
Problem 1.
Write the corresponding DE.
In an interval ∆t the loss and gains are
∆Q = ∆t { -R + K ( Cext Cint (t) ) }
Dividing by (∆t V) wher v is the cells volume we get
dC/dt = { -R/V +(K/V) ( Cext Cint (t)) }.
When the steady state value is achieved dC/dt=0 ,so
Cext Cint = (V/K)(R/V) or
Cint = Cext R/K .
Let -R/V +(K/V)Cext = a and (K/V) = b . The Laplace transform with C(0)=c0 of the DE gives
sc(s) c0 =a/s bc(s) , thus ;
c(s) = c0/(s+b) +a/(s*(s+b))
Employing MATLAB we get
syms a b s t c0 ;
>> cs=c0/(s+b) +a/(s*(s+b)) ;
>> ilaplace(cs)
ans =
C(t) = a/b+exp(-b*t)*(c0*b-a)/b .
Problem 4. An organ and its contained blood are represented by a two compartment scheme , see figure. Both contain identical initial concentrations of a tracer C0 = C1(0) = C2(0). There is diffusion of the tracer between blood and tissue. At time zero blood begins to flow outside of tissue. Find the DE for C1 and C2 and show that C2 satisfies the 2nd order
V1 V2 d2 C2/dt2 +[V1K +(F+K)V2 ] dC2/dt + FKC2 = 0 (a)
Figure for problem 4.
Eq (a) is corresponds to a damped oscillator ,see for example eq. 9 table III-1 of Simon textbook. Let in this case K12 < 4 K2. Part of the exponential solution become imaginary i.e. produces oscillatory motion.
Our solution:
Recall always that C~moles/volume , K ~volume/time
F~ volume/time
From what was explained in previous lectures about the two compartment problem we pose the DE ,
V1 d C1/dt = - K (C1 C2) F C1 ~ moles/time (b)
V2 d C2/dt = +K (C1 C2) ~ moles/time (c )
with IC C0 = C1(0) = C2(0).
In the next plot we have solved numerically (b) and (c) simultaneously with arbitray values C1(0)=C2(0)= 0.5 , V1=V2=1 , K=0.1 ,F=.05 . They do not correspond to any particular physiological problem.
Fig . Solutions of eqs (b) (c).
FORTRAN code solution (b) (c)
c problem 4 page 110 Math Tech for Biology and Medicine
c V1 d C1/dt = - K (C1 - C2) - F C1
c V2 d C2/dt = +K (C1 - C2)
real k
data c1zero, c2zero /.5 ,.5 /
data V1,V2 , K , f / 1. ,1., .1 ,.05/
a1(c1,c2)=(1./v1)*( -k*(c1-c2) -f*c1)
a2(c1,c2)=(1./v2)*(k*(c1-c2))
tsmall= amin1(V1/k,V2/k,V1/f)
tlarge=amax1(V1/k,V2/k,V1/f)
dt=tsmall/100.
tf=5.*tlarge
nstep=int(tf/dt)
kp=int(float(nstep)/60.)
kount=kp
print 100 ,0. , c1zero,c2zero
do 10 i=1,nstep
t=dt*float(i)
c1one=c1zero + dt*a1(c1zero,c2zero)
c2one=c2zero+ dt*a2(c1zero,c2zero)
if(i.eq.kount)then
print 100 , t, c1one,c2one
kount=kount+kp
endif
c1zero=c1one
c2zero=c2one
10 continue
100 format(1x,'t,c1,c2=',3(3x,e10.3))
stop
end
The numerical solution of eq(a) , the second order eq for C2 is worked out next.
IC , C2(0)=0.5 and from (b) dC2(0)/dt=0.
c problem 4 page 110 Math Tech for Biology and Medicine
c V1 d C1/dt = - K (C1 - C2) - F C1
c V2 d C2/dt = +K (C1 - C2)
c V1 V2 d^2c2/dt^2 = -[V1*K +(F+K)*V2 ]* dC2/dt - F*K*C2
real k
data c1zero, c2zero /.5 ,.5 /
data V1,V2 , K , f / 1. ,1., .1 ,.05/
a3( c2zero,c2one)=(-1./(v1*v2))*( (v1*k+(f+k))*(c2one-c2zero)/dt
$ + f*k*c2one)
dc2dt=(1./v2)*k*(c1zero-c2zero)
tsmall= amin1(V1/k,V2/k,V1/f)
tlarge=amax1(V1/k,V2/k,V1/f)
dt=tsmall/100.
c2one=c2zero+dt*dc2dt
tf=5.*tlarge
nstep=int(tf/dt)
kp=int(float(nstep)/60.)
kount=kp
print 100 ,0. , c2zero
do 10 i=1,nstep
t=dt*float(i)
c2two=2.*c2one-c2zero+dt**2*a3(c2zero,c2one)
if(i.eq.kount)then
print 100 , t, c2two
kount=kount+kp
endif
c2zero=c2one
c2one=c2two
10 continue
100 format(1x,'t,c2=',2(3x,e10.3))
stop
end
Fig . Solution C2 from the 2nd order DE , eq (a).
We show now how to go from the first order set (b)-(c) to (a).
Take the time derivative of V1 * eq(c)
V1 V2 d2 C2/dt2 = +K V1 ( dC1/dt dC2 /dt )
= K V1 dC1/dt
- K V1 dC2/dt . (d)
Using eq (b) in (d)
V1 V2 d2 C2/dt2 = K{- K (C1 C2) F C1 }
- K V1 dC2/dt
= { -KV2 dC2/dt KF C1} - K V1 dC2/dt
= - { KV1 + KV2 } dC2/dt -FKC1
Finally substitute for KC1 = V2 dC2/dt + KC2 to arrive at eq (a),
V1 V2 d2 C2/dt2 = - { KV1 + KV2 } dC2/dt -F {V2 dC2/dt + KC2}
= -{ FV2 + K(V1 + V2) } dC2/dt - FKC2
This shows once that it is possible to go from a set of two 1rst order DE to a second order DE.
However we showed before that we can deal numerically with the simultaneous set.
END OF LECTURE 5