Lecture 5.

 

Compartmental Problems

 

by Reinaldo Baretti Machn

 

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.

 

 

Section 2. The one compartment mixed case

 

Suppose that into one compartment I ~moles/time are flowing and at the same time there is a loss , a diffusion of L ~ volume/time.

A quantity of moles Q is changing in an interval ∆t by

 

∆Q = ∆t I   - ∆t  L C(t)                                         (1)

 

moles= time*(moles/time) time*(volume/time)*(moles/volume).

 

The first term , ~I  is independent of the concentration. That is why it is called a mixed case. With appropriate value of  I , C and  L  this may represent a cell from which a susbtance diffuses out while theres is a constant input I.

 

 

 

Dividing by V and taking the limit as ∆t goes to zero gives the DE

 

dC(t)/dt  = I/V (L/V) C(t)                                . (2)

Thre cases are posed.

First the steady state, when dC/dt becomes zero, has the solution ,

 

Csteady state = I/V                                                 (3)

 

Second suppose the initial concentration C(0)=0. Solve (2) using the Laplace transform. From the previous chapters we anticipate that the answer is an inverted exponential.

 

From

L { dC(t)/dt  = I/V (L/V) C(t)  }

sc(s) C(0) = I/(s*V) (L/V) c(s)     ,

we find

 

c(s) = (I/V)* (1/(s*(s+L/V)))  

 

Using Matlab

>> syms I V L s t;

>> cs=(I/V)* (1/(s*(s+L/V)));

>> ilaplace(cs)

 

ans =

 

I/L*(1-exp(-L*t/V))

 

So C(t) = (I/L) { 1 exp(-Lt/V) }

 

Let  I=.01 moles/minute  , L=0.2 liters/minute ,  C(0)=0.03 moles/liter,

V=10 liters .

 

I=.01 ; L=0.2 ; c0=.3 ;V=10 ;

tscale=V/L ;tf=3*tscale ;

t=[0:tf/200:tf];

ct=(I/L)*( 1 - exp(-L*t/V) );

>> plot(t,ct),xlabel('t~minutes'), ylabel('c~moles/liter')

>> 

Fig 5.1 . An inverted exponential results for C(t).

 

 

 

Lecture 3. The Two-Compartment Problem

 

Fig 5.2 The two compartment problem.

 

The meaning and dimensions of the terms are ,

C1(t) =Q1/V1   , C2(t) = Q2/V2 ~moles/volume

 

I1 and I2 ~ moles/time are inputs to the compartments

 

K(C1-C2) ~ diffusion between compartments

            ~  (volume/time)*(moles/volume)

 

L1C1 , L2 C2 ~ concentration dependent loss

                    ~ (volume/time)*(moles/volume)          

In an interval ∆t the net changes in Q1 and Q2 are ,

 

∆ Q1 = ∆t { I1 - L1C1 - K(C1-C2) }                                    (4)

 

∆ Q2 = ∆t { I2 L2C2 + K(C1-C2) }                                   (5)

 

Notice that the last term in (5) must have the opposite sign of the same term appearing in(4).

Dividing (4) by ∆t V1 and (5) by ∆t V2 give the two DE,

 

dC1/dt = (1/V1)*{ I1 - L1C1 - K(C1-C2) }            (6)

 

d C2/dt = (1/V2)* { I2 L2C2 + K(C1-C2) }          (7)

 

 

We solve (6) & (7) using finite differences with the following data and initial conditions (IC). This is not an example from physiology. The values of the parameters are chosen to illustrate a solution.

L ~liters/minute , K~liters/minute , I ~ moles/minute , C~moles/liter,

V~liters

I1=.01;I2=.02 ;L1=.1 ;L2=.2;Kd=.5*L1;

c1(1)=0.20 ; c2(1)=0.1 ;

V1=10 ; V2= 10 ;

 

 

 

Fig 5.3 Concentrations C1 , C2 . C1 is the upper curve.

 

% chapter5 sec3

I1=.01;I2=.02 ;L1=.1 ;L2=.2;Kd=.5*L1;

c1(1)=0.20 ; c2(1)=0.1 ;

V1=10 ; V2= 10 ;

A=[V1/L1 V1/Kd V2/L2 V2/Kd];

tscale=min(A);tlarge=max(A);

tf=tlarge; dt=tscale/50;

k=0;

for t=[dt:dt:tf];

    k=k+1;

    g1=(I1-L1*c1(k)-Kd*(c1(k)-c2(k)))/V1;

    g2=(I2-L2*c2(k)-Kd*(c2(k)-c1(k)))/V2;

    c1(k+1)=c1(k)+dt*g1;

    c2(k+1)=c2(k)+dt*g2;

end 

t=[0:dt:tf];

plot(t,c1,t,c2),xlabel('t~minutes'),...

    ylabel('c1,c2~moles/liter')

 

 

We work next a case  proposed in page 89 of Reference  1.

 

Consider the case I2 =0 , C1(0)=C2(0)=0.

MATLAB code

% chapter5 sec3

I1=.01; I2=.0 ;L1=.1 ;L2=.2; Kd=.5*L1;

c1(1)=0.0 ; c2(1)=0.0 ;

V1=10 ; V2= 10 ;

A=[V1/L1 V1/Kd V2/L2 V2/Kd];

tscale=min(A);tlarge=max(A);

tf=tlarge; dt=tscale/50;

k=0;

for t=[dt:dt:tf];

    k=k+1;

    g1=(I1-L1*c1(k)-Kd*(c1(k)-c2(k)))/V1;

    g2=(I2-L2*c2(k)-Kd*(c2(k)-c1(k)))/V2;

    c1(k+1)=c1(k)+dt*g1;

    c2(k+1)=c2(k)+dt*g2;

end 

t=[0:dt:tf];

plot(t,c1,t,c2),xlabel('t~minutes'),...

    ylabel('moles/liter')

                                               

    *     

Fig. 5.4 Concentrations C1 , C2 .   

 

 

 

 

 

 We work now  a second  case  proposed in page 91 of Reference  1.

Let   K=3*(L1 + L2) and   L2 << L1 . The system will behave as if it were a single compartment.

 

Data   I1=.1  , I2 =0 , C1(0)=C2(0)=0.

 

L1=.2 ;L2 =.05=.2 ; K=3*(L1+L2)  , ,V1 =V2 =5

 

 The units are again, L ~liters/minute , K~liters/minute , I ~ moles/minute , C~moles/liter,  V~liters.

 

 

Fig 5.5 Concentration of C1 , C2 with a large diffusion from the first compartment  to the second.  

 

  

 

MATLAB code

 I1=.1; I2=.0 ; L1=.2 ; L2=.05 ;

Kd=3*(L1+L2); c1(1)=0.0 ; c2(1)=0.0 ;

V1=5 ; V2= 5 ;

A=[V1/L1 V1/Kd V2/L2 V2/Kd];

tscale=min(A);tlarge=max(A);

tf=tlarge; dt=tscale/50;

k=0;

for t=[dt:dt:tf];

    k=k+1;

    g1=(I1-L1*c1(k)-Kd*(c1(k)-c2(k)))/V1;

    g2=(I2-L2*c2(k)-Kd*(c2(k)-c1(k)))/V2;

    c1(k+1)=c1(k)+dt*g1;

    c2(k+1)=c2(k)+dt*g2;

end 

t=[0:dt:tf];

plot(t,c1,t,c2),xlabel('t~minutes'),...

    ylabel('moles/liter')      

 

 

Section 4. A simple three compartment problem     

Fig 5.6. Three compartments with a constant input I3 at #3 and with outflow

L C3(t) .There is diffsusion between compartments.

 

 

Let C~moles/minute    , V ~liters  , L~ liters/minute , K~liters/minute ,

and I ~ moles/minute.

 

The system is described by three coupled first order linear DE.

 

V1 dC1/dt = -K12 ( C1 C2)                                                         (8)

 

liter*(moles/(liter-miniute)) = (liter/minute)*(moles/liter)

 

Compartment #2 has an incoming diffusion  + K12 ( C1 C2) and the

 

outgoing diffusion   -K23 (C2 C3)   . The DE is                                                    

 

V2 dC2/dt = + K12 ( C1 C2)   -K23 (C2 C3)            .     (9)   

 

The parameters  K12  , K23    , need not be the same. 

Compartment #3 has a steady input I3  , an ingoing diffusion

 

 +K23 (C2 C3)   and the out flow LC3 . The DE is

 

V3 dC3/dt = I3  + K23 (C2 C3)   -L C3                    .   (10) 

 

The following data is used for the next plot.

  Data

I3=.2; L3=.05 ; ;K12=.01 ;K23=.01;

c1(1)=0.50 ; c2(1)=0.;c3(1)=0. ;

V1=5 ; V2= 5 ;V3=5;

 

 

 

 MATLAB

 % chapter5 sec4 Three compartment

I3=.2; L3=.05 ; ;K12=.01 ;K23=.01;

c1(1)=0.50 ; c2(1)=0.;c3(1)=0. ;

V1=5 ; V2= 5 ;V3=5;

A=[V1/L3 V1/K12 V1/K23 V2/K23 V2/L3 V3/L3];

tscale=min(A);tlarge=max(A);

tf=tlarge; dt=tf/3000;

k=0;

for t=[dt:dt:tf];

    k=k+1;

    g1=-K12*(c1(k)-c2(k))/V1;

    g2=(K12*(c1(k)-c2(k))-K23*(c2(k)-c3(k)))/V2;                                                   

    g3=(I3-L3*c3(k)+K23*(c2(k)-c3(k) ) )/V3;

    c1(k+1)=c1(k)+dt*g1;

    c2(k+1)=c2(k)+dt*g2;

    c3(k+1)=c3(k)+dt*g3;

end 

t=[0:dt:tf];

plot(t,c1,t,c2,t,c3),xlabel('t~minutes'),...

    ylabel('c1,c2,c3~moles/liter')

 

 

 

 

       

 

                                   

       

Fig 5.6 Concentrations C1 , C2 , C3 in the three compartment problem.  

 

 

 

Section 5.  The delta impulse function

 

One  method for solving inhomogeneous  DE is that of the Greens function.

 

We explain the method through a concrete example

 

Let the DE be

 

d2 C/dt2 +  dC/dt +  C = F(t)            , t ≥ 0              (11)

 

with F(t) = cos(2t)    and IC  C(0)=1 , dC(0)/dt=0.

 

The complete solution of  (11) has two parts , the particular solution Cp(t)

and the complementary or homogeneous solution  Ch(t) .

 

C(t)  = Cp(t) + Ch(t)                                             . (12)

 

The homogenous part satisfies

 

d2 Ch /dt2 +  dCh /dt +  Ch = 0   and the particular solution satisfies

 

   d2 Cp /dt2 +  dCp/dt +  Cp = F(t)    .

 

The Green function satisfies the DE with a delta function on the right hand side             

 

 d2 G /dt2 +  dG/dt +  G= δ (t)                                          .   (13)

 

δ (t)   is  Dirac delta function .Some of its properties are

 

∫  δ (t)  dt = 1   ,       - ∞   <t < ∞      ,

 

∫ F(t)  δ (t-a)  dt    = F(a)  ,    - ∞   <t < ∞        .                                                              

 

Once G(t) is obtained , the particular solution is given by the  integral

 

Cp(t) = ∫t 0 G(t-t) F(t) dt

 

The complete solution can be written as

 

                       C(t) = ∫t 0 G(t-t) F(t) dt  + Ch(t)   .                 (14)

 

Taking the derivatives upon (14) we verify that

 

d2 C(t) /dt2 +  dC(t)/dt +  C(t) =

 

∫ {d2 G(t-t) /dt2 +  dG(t-t)/dt +  G(t-t)} F(t) dt

 

+  d2 Cp(t) /dt2 +  dCp(t)/dt +  Ch(t)

 

= ∫   δ(t-t) F(t) dt   + 0 =  F(t)                                             . (15)

 

 

Example :

We solve (11) in stages. First find the homogeneous solution.

 

Part 1. Solution to the homogenous equation

d2 C/dt2 +  dC/dt +  C =0  with IC , C(0)=1, dC(0)/dt=0 is worked using the

Laplace transform.

 

L (d2 C/dt2 +  dC/dt +  C) = s2 c(s) sC(0)-dC(0)/dt +sc(s)-C(0) +c(s)=0 .

 

c(s) = (s+1)/(s2 +s +1)

MATLAB code

>> syms s t;

>> cs=(1+s)/(s^2+s+1);

>> ilaplace(cs)

 

gives

 

C(t) = exp(-1/2*t)*cos(1/2*3^(1/2)*t)+1/3*exp(-1/2*t)*3^(1/2)*sin(1/2*3^(1/2)*t)

 

 

 

t=[0:.05:10];

% Gcomp is the complementary solution with IC g(0)=1

% dG(0)/dt=0

Gcomp=exp(-t./2).*cos(sqrt(3)*t./2) + ...

sqrt(3)*(1/3)*exp(-t./2).*sin(sqrt(3)*t./2);

plot(t,Gcomp),xlabel('t'),ylabel('Gcomplementary')

 

Fig 5.7    . Complementary solution C comp. ( or Chomogeneous ) , i.e. the solution of the homogeneous DE.

 

Part 2.  Solution to the delta function impulse

d2 G/dt2 +  dG/dt +  G = δ(t)     , with G(0)=0 , dG(0)/dt=0

 

s2 g(s) sG(0)-dG(0)/dt + sg(s) G(0) + g(s) = 1       ( 16  )

 

g(s) = 1/(s2 + s + 1)                                            ( 17  )

 

MATLAB

syms s t;

>> f=1/(s^2+s+1);

>> ilaplace(f)

 

ans = 2/3*3^(1/2)*exp(-1/2*t)*sin(1/2*3^(1/2)*t)

 

G(t)=(2/31/2)exp(-t/2) sin(31/2 t/2)

and

G(t-t) = (2/31/2)exp(-(t-t)/2) sin(31/2 (t-t)/2)

 

We plot the solution to the delta function impulse , the Green function for the DE.

MATLAB

t=[0:.05:10];

G=(2/3^(1/2))*exp(-t./2).* sin(3^(1/2)* t./2);

plot(t,G),xlabel('t'),ylabel('G')

 

 

 

 

Fig   5.8   . Greens function.

 

 

As was said before, the solution to

 

d2 C/dt2 +  dC/dt +  C = F(t)             t ≥ 0 is

 

 

C(t) = ∫ t 0 G(t-t) F(t) dt  + C(t)complementary      

The following FORTRAN code calculates the integral  

t 0 G(t-t) F(t) dt  and adds  C(t)complementary  . The plot is given .

  

 

 

Fig 5.9  Plot of C(t) = ∫ t 0 G(t-t) F(t) dt  + C(t)complementary    

 

FORTRAN code

c Use of Green's function to solve inhomogenous  DE 22 nov 2008

      data ti ,tf ,nt, nstep/0.,10.,70, 5000/

c g(u=t-t') is the greens'function with L g =delta(t) g(0)=0,dg/dt=0.

      g(u)=(2./3.)*sqrt(3.)*exp(-u/2.)*sin(sqrt(3.)*u/2.)

c complentary solutio to L C = 0 , C(0)=1. , dC/dt=0.

      comp(t)=exp(-t/2.)*cos(sqrt(3.)*t/2.)+

     $ (1./3.)*exp(-t/2.)*sqrt(3.)*sin(sqrt(3.)*t/2.)

c f(t)  is the driving force

      f(t)=cos(2.*t)

      dt=(tf-ti)/float(nt)

      do 10 it=0,nt

      t=ti+dt*float(it)

      dtprime=t/float(nstep)

      du=dtprime

      sum=0.

      do 20 j=1,nstep

      tprime=dtprime*float(j)

      u=t-tprime

      sum=sum+(dtprime/2.)*(g(u)*f(tprime)+g(u-du)*f(tprime-du))

20    continue

      ctotal=sum+comp(t)

      print 100 ,t,sum,comp(t),ctotal

10    continue

100   format(1x,'t,cinh,comp,c(t)=',4(3x,e10.3))

      stop

      end

 

We end this section by solving the problem by the finite difference method.

The generic solution to d2 C/dt2 +  dC/dt +  C = F(t)           is simply

 

C (tn ) = 2 C(tn-1) C(tn-2) + (∆t)2 { - (C(tn-1)-C(tn))/ ∆t C(tn-1) + F(tn-1)  }

The  IC is C(0)=1 , dC(0)/dt=0 which is equivalent to   C(t0) =1 , C(t1) =1  . 

 

 

The solution is plotted.

Fig 5.10 . C(t) obtained by finite differences. It is identical to fig 5.9.

 

Fortran code  Solution by finite differences

c23 nov 2008 equation d^2c/dt^2+dc/dt+c= cos(2t) IC c(0)=1 ,dc(0)/dt=0.

c solution by finite differences

      data ti, tf , nstep/0.,10.,4000 /

      data c0,c1/1.,1./

      f(t)= cos(2.*t)

      a(t)=-(c1-c0)/dt -c1 +f(t)

      dt=(tf-ti)/float(nstep)

      kp=int(float(nstep)/70.)

      kount=kp

      print 100,0.,c0

      do 10 i=2,nstep

      t=ti+dt*float(i)

      c2=2.*c1-c0+dt**2*a(t-dt)

      if(i.eq.kount)then

      print 100,t,c2

      kount=kount+kp

      endif

      c0=c1

      c1=c2

10    continue

      print 100,t,c2

100   format(1x,'t,C(t)=',2(4x,e10.3))

      stop

      end