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:

reibaretti2004@yahoo.com

 

 

 

 

 

 

Reference :

 

1. Mathematical Techniques for Biology and Medicine (Dover Publications- Paperback)

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

 

7. http://www.mathworks.com/access/helpdesk/help/toolbox/symbolic/index.html?/access/helpdesk/help/toolbox/symbolic/laplace.html&http://search.yahoo.com/search?p=MATLAB+laplace+transfom&ei=utf-8&fr=b1ie7

 

 

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 cell’s 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