Lecture 7.  

 

Regulation and oscillation

 

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 7 of  Reference 1. The problem numbers refer to the textbook listing.

 

Section 3. A simple regulator.

 

Fig 1. By means of a float the volume level  V is monitored. The feedback mechanism is sketched. The rate Ra is accordingly increased or decreased.

 

 As simple regulator is shown in figure 1.

 

The rate of change of V can be expressed with a first order DE.

Let VZF be the reference volume. If V= VZF the feedback mechanism cut Ra to zero. If V < VZF pump A should open the flow Ra.

 

 

Then

              dV/dt = G (VZF –V) – Rb                      (1)

where G ~1/time      and   Rb ~ volume/time    . We assume Rb is constant.

 

Eq (1) is of the type already encountered

 

dy/dy = a – b y          with IC y(0).

Let y(t) = V(t) , a  = G VZF – Rb   and b = G .The solution is

 

y(t) = a/b + ( y(0) –a/b) exp(-bt)  ,   or  

 

V(t) = ( VZF – Rb /G ) + {V(0) – (  VZF – Rb /G ) } exp(-Gt)   (2)

 

                  

 

 

Example: (These are not typical physiological values). Let  V(0)= 12 liters,

VZF =10  liters , G = 1/(2 minutes)  = .5 min-1    Rb = 1.0 liter/min

 

The finite difference solution is

 

V(tn) = V(tn-1) + ∆t {  G (VZF –V(tn) ) – Rb  }         (3)

 

The time scale is 1/G  = 2 minutes.

 

The steady state solution is achieved when

   G (VZF –V(tn) ) – Rb    = 0 ,

or    Vsteady state = VZF  -Rb/G = 10-1/.5=8 liters ( see next figure)

 

Fig 2 .Volume with simple regulator mechanism.

 

MATALAB code

% Simon page 135

g=.5 ; v(1)=12 ;vzf=10 ;Rb=1;

tscale=1/g; dt=tscale/200 ;tfinal=4*tscale;

k=1;

for t=[dt:dt:tfinal];

    k=k+1;

    a=g*(vzf-v(k-1)) -Rb;

    v(k)=v(k-1)+dt*a;

end

 

t=[0:dt:tfinal];

plot(t,v),xlabel('t~minutes'),ylabel('V~liters')

 

 

 

 

Suppose now we have V(0)= Vsteady state = 8 liters and at t=2 minutes a quantity Vimpulse = 2 liters is dropped in the compartment. Write the corresponding DE and solve. Equation (1) is modified with the introduction of the delta function

 

              dV/dt = G (VZF –V) – Rb   + 2 δ (t-2)              (4)

 

 

          We use again the finite difference solution

 

   V(tn) = V(tn-1) + ∆t {  G (VZF –V(tn) ) – Rb  }     2 δ (tn-2)    (5)

 

We define

 

  δ (tn-2)    = 1/(∆t)     ,              2  ≤  tn  ≤  2 + ∆t    ,

 

                 =   0                           for any other t value.  

 

 

FuiFig

Fig 3. Tank volume with impulse disturbance at t=2 minutes.

 

FORTRAN code

c W. Simon page 137

      data g,v0,vzf,rb/.5,8.,10.,1./

      tscale=1./g

      dt=tscale/200.

      tfinal=8.

      nstep=int(tfinal/dt)

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

      kount=kp

      print 100 ,0.,v0

      do 10 i=1,nstep

      t=dt*float(i)

      a=g*(vzf-v0) -rb +2.*delta(t-2.,dt)

      v1=v0+dt*a

      if(i.eq.kount)then

      print 100 ,t,v1

      kount=kount+kp

      endif

      v0=v1

10    continue

100   format(1x,',t(min),v(L)=',2(3x,e10.3) )

      stop

      end

     

      function delta(u,dt)

      if(u.ge.0..and.u.le.dt)delta=1./dt

      if(u.ne.0.)delta=0.

      return

      end

 

 

 

Section 5. The effect of time delay in the feedback loop

 

Suppose the height of pump A , in Fig. 1, is increased such that the time of fall tdelay becomes an appreciable magnitude of the order of 1/G. Then the pump A is responding to a volume value at a previous instant equal to t-tdelay .

 

Equation (1)  dV/dt = G {VZF –V(t) } – Rb  , should be modified to

 

take into account such delay  ,

 

                  

                              dV/dt = G {VZF –V(t-tdelay) } – Rb  .            (6)

 

Fig 4. The effect of a time delay of 1 sec when G= 2.1 /sec . The amplitude of the oscillation is unbounded.

 

 

 

Fig 5. The effect of a time delay =1/G=  0.476 sec . The amplitude of the oscillation is bounded and attenuated.

 

 

FORTRAN code for eq (6)

c W. Simon page 137

      dimension v(0:10000)

      data g,vzf,rb/2.1,10.,1./

      tscale=1./g

      tdelay=1.

      tsmall=amin1(tdelay,tscale)

      dt=tsmall/200.

      tfinal=15.

      nstep=int(tfinal/dt)

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

      kount=kp

      v(0)=12.

      id=int(tdelay/dt)

      print 100 ,0.,v(0)

      do 10 i=1,nstep

      t=dt*float(i)

      if(i.le.id)vdelay=v(0)

      if(i.gt.id)vdelay=v(i-id)

      a=g*(vzf-vdelay) -rb

      v(i)=v(i-1)+dt*a

      if(i.eq.kount)then

      print 100 ,t,v(i)

      kount=kount+kp

      endif

10    continue

100   format(1x,',t(min),v(L)=',2(3x,e10.3) )

      stop

      end

     

 

Section 6. Pooling delay

 

 

 

 

 

Fig 6. Pump A receives feedback signal from Volume 2. V2 is filled not directly by RA but in a delayed manner at the rate K1V1 .

 

The coupled DE’s  are

 

                              dV1/dt = RA – K1 V1                      (7)

but   RA is now a variable controlled by the feedback , thus

              RA(t)  = G ( VZF – V2 (t) )     and

 

      dV1/dt = G(VZF – V2 )  – K1 V1                           (8)

 

with G~1/time ,  K1 ~ 1/time.

 

For the second volume

 

dV2/dt = K1 V1 – RB                                                (9)

 

where  RB ~ volume/time and it is assumed to be constant.                                                    

 

The difference equations are

 

      a1=g*(vzf-v20)-k1*v10

      a2=k1*v10-rb

      v11=v10+dt*a1

      v21=v20+dt*a2

or

 

V1 (tn) = V1 (tn-1)  + ∆t {  g(VZF – V2 (tn-1 ) –K1 V1(tn-1) }

 

 

 

Example: Fig 7 shows the oscillations of the volumes with using   the data g=4 /s , vzf =10 cm3 , k1=.1/s , rb=1cm3/s  , v1(0)=8 cm3 ,v2(0)= 11 cm3 .

Figure 7. Volume in the two tanks with pooling delay. V 2 is the purple curve with initial value of 11.

 

 

 

 

Example 2 (fig 8 )    data g ,vzf,k1/.2,10.,.05/

      data rb,v10,v20/1.,30.,8./

Fig 8. Oscillations with second set of data.

 

 

FORTRAN code for figure 7.

c Simon page 144

      real k1

      data g ,vzf,k1/4.,10.,.1/

      data rb,v10,v20/1.,8.,11./

      data ti/0./

      tscale=amin1(1./g,1./k1)

      tfinal=5.*amax1(1./g,1./k1)

      dt=tscale/200.

      nstep=int(tfinal/dt)

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

      kount=kp

      print 100,0., v10,v20

      do 10 i=1,nstep

      t=ti+dt*float(i)

      a1=g*(vzf-v20)-k1*v10

      a2=k1*v10-rb

      v11=v10+dt*a1

      v21=v20+dt*a2

      if(i.eq.kount)then

      print 100,t, v11,v21

      kount=kount+kp

      endif

      v10=v11

      v20=v21

10    continue

100   format(1x,'t,v1,v2=',3(3x,e10.3))

      stop

      end

     

 

 

An arbitrary set of data for G  , K1 , RB , VZF , V1(0) and V2(0) may easily produce negative values of V1 and V2 which may be interpreted as a non physical situation.

 

 

 Examples 4 & 5. An important problem of physiological importance are the oscillations of the insulin storage in the pancreas and storage in the blood.

Mathematically the problem would be somewhat like what happens in these two examples.

 

 

In figure 9 a set of values produces overshoot or oscillations above and below the steady state values.

 

In Figure 10 the oscillations are eliminated by decreasing RB from 20 cm3/s    to 12 cm3/s   .  decreasing G from 0.25/s to .025/s and increasing VZF from

240 to 600 cm3. The decrease in G is equivalent to a larger delay time.

 

 

 

Fig 9.  A given set of parameters produces oscillation or overshoot above and below the steady state values.

 

 

Fig 10. The various parameters have been adjusted to eliminate oscillations.

 

 

FORTRAN code

c Simon page 150 topic of overshhot

      real k1

      data g ,rb,k1/.25,20.,.1/

      data vzf,v10,v20/240.,100.,200./

      data ti/0./

c      vzf=k1*200./g +160.

      print*,'vzf=',vzf

      print*,'  '

      tscale=amin1(1./g,1./k1)

      tfinal=6.*amax1(1./g,1./k1)

      dt=tscale/250.

      nstep=int(tfinal/dt)

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

      kount=kp

      print 100,0.,g*(vzf-v20), v10,v20

      do 10 i=1,nstep

      t=ti+dt*float(i)

      a1=g*(vzf-v20)-k1*v10

      a2=k1*v10-rb

      v11=v10+dt*a1

      v21=v20+dt*a2

      ra=g*(vzf-v20)

      if(i.eq.kount)then

      print 100,t, ra,v11,v21

      kount=kount+kp

      endif

      v10=v11

c      if(v10.lt.0.)v10=0.

      v20=v21

c      if(v20.lt.0.)v20=0.

10    continue

100   format(1x,'t,ra,v1,v2=',4(3x,e10.3))

      stop

      end

     

     

  

Another way of controlling the overshoot is to include  a term –M dV2/dt in

the definition of  Ra. This is called a velocity feedback , and we have

 

               RA(t)  = G ( VZF – V2 (t) )   - M dV2/dt  .

 

The FORTRAN code  ( see code below)is modified to include M and redefine

    

      a2=k1*v10-rb

      a1=g*(vzf-v20)-M*a2-k1*v10

 

 

 

Figures 11 and 12 show the volumes oscillations for two different

 values of M .

 

Fig 11. Damped flow with velocity feedback (M=0.1).

 

 

 

Fig 12. The overshoot of Fig 11 is reduced by increasing the

                   parameter M to 1.0 .

 

FORTRAN code

c FORTRAN code for Simon page 150 introduction of -M dV2/dt term in Ra

      real k1 ,M

      data g ,m,vzf,k1/.25,.1,240.,.1/

      data rb,v10,v20/20.,100.,200./

      data ti/0./

      tscale=amin1(1./g,1./k1)

      tfinal=5.*amax1(1./g,1./k1)

      dt=tscale/200.

      nstep=int(tfinal/dt)

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

      kount=kp

      print 100,0., v10,v20

      do 10 i=1,nstep

      t=ti+dt*float(i)

c      a1=g*(vzf-v20)-M*a2-k1*v10

      a2=k1*v10-rb

      a1=g*(vzf-v20)-M*a2-k1*v10

      v11=v10+dt*a1

      v21=v20+dt*a2

      if(i.eq.kount)then

      print 100,t, v11,v21

      kount=kount+kp

      endif

      v10=v11

      v20=v21

10    continue

100   format(1x,'t,v1,v2=',3(3x,e10.3))

      stop

      end

 

 

 

Section 7.  A Three Compartment System

A three compartment system presents many possible options , one of which is shown in Figure 13. But one could have for example the level  V3 serving a feedback to pump A etc..

 

 

Fig 13. A three compartment system in which the level V3 controls the flow

 

RC . The parameter RA would be taken as a function of time while RB would be assumed constant.

 

 

 

Let RA vary armonically with frequency ω ,     

   

   RA (t)= RA0  + A sin( ω t)  ~ cm3 /s                     .                   (7.1)

 

RC   is controlled by the feedback from V3 ,

                     RC (t)  = G ( VZF – V3(t) )  ~ cm3 /s                      (7.2)

where   G ~ 1 /seconds .

 

The three DE , one for each volume are , (all Ki ~1/sec )

 

                    dV1/dt = RA – K1 V1                ,             (7.3)

 

                    dV2/dt = K1 V1 – K2 V2          ,             (7.4)

 

 

                   dV3/dt = K2 V2 – RB                 ,            (7.5)        

 

 One must have the initial conditions ,   V1 (0)  ,  V2 (0),  V3 (0) .                                                                     

Three coupled DE would lead ,after some substitutions, to third order DE for each volume , see Ref 1 page 153. 

But we can work 7.3 – 7.5 by finite differences. The solutions are

 

 

V1(tn)  = V1(tn-1)  + (∆t) {RA (tn-1)  – K1 V1 (tn-1)  }  ,      (7.6)

 

 V2(tn)  = V2(tn-1)  + (∆t) { K1 V1 (tn-1) – K2 V2 (tn-1) }  ,  (7.7) 

 

V3(tn)  = V3(tn-1)  + (∆t) { K2 V2 (tn-1) – RB  }       .          (7.8)        

 

 

There are for scales of time in this problem and the outcome can depend on their relative magnitudes.

 

One has  tG = 1/G , t K1 = 1/K1 , tK2 = 1/K2  , tω = 2π/ω . The interval

 

∆t has to be much smaller then the smallest of the time scales.

 

 

  

Fig 14 . Three compartment problem.

 

Fig 15.  RC 

c FORTRAN code for Simon page 152 three compartments

      real k1 ,k2

      data g ,vzf,k1,k2/.25,260.,.1,.1/

      data rb,v10,v20,v30/20.,100.,200.,250./

      data ti/0./

      ra(t)=20.+ 5.*sin(w*t)

      pi=2.*asin(1.)

      w=2.*pi*G

      tscale=amin1(1./g,1./k1,1./k2,2.*pi/w)

      tfinal=5.*amax1(1./g,1./k1,1./k2,2.*pi/w)

      dt=tscale/200.

      nstep=int(tfinal/dt)

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

      kount=kp

      print 100,0., v10,v20,v30,g*(vzf-v30)

      do 10 i=1,nstep

      t=ti+dt*float(i)

      a1= ra(t-dt)-k1*v10

      a2=k1*v10-k2*v20

      a3=K2*V20 - RB

      v11=v10+dt*a1

      v21=v20+dt*a2

      v31=v30+dt*a3

      rc=g*(vzf-v31)

      if(i.eq.kount)then

      print 100,t,v11,v21,v31,rc

      kount=kount+kp

      endif

      v10=v11

      v20=v21

      v30=v31

10    continue

100   format('t,v1,v2,v3,rc=',5(3x,e10.3))

      stop

      end

 

 

 

 

Section 8. Control of Biosynthesis

 

A biological chemical reaction involving an enzyme is described by the diagram

 

(8.1)

Fig 16. Reaction  involving an enzyme.

The meaning of the symbols an their units are

S substate ~ moles/volume

 

EF ~ free enzyme ~ moles/volume

 I~ intermediate product ~ moles/volume

P ~  product ~ moles/volume

 

The free enzyme is related to the total enzyme ET by

 

                EF (t)  = ET – I(t)                                                  (8.2)

 

The kinetics of the reaction are represented by

 

d I/dt = k1 S EF – (k2+k3) I                                     ,           (8.3)

 

or

 

d I/dt = k1 S (ET –I ) – (k2+k3) I                      ,                  (8.4)                           

 

where k1~ (volume/(mol-time))    ,   k2 ~ 1/time   , k3 ~ 1/time.

 

d P /dt = k3 I                                          .                          (8.5)

 

The initial conditions are assumed to be I(0)=P(0)=0.

To solve (8.4) –(8.5) we supposes that k1 ,k2 k3 are given ands S and ET are constants.

  The Michaelis constant  KM ~ moles/volume is defined by

                

                             KM ≡ (k2+k3)/k1                            (8.6)

 

The substrate S is large or small when compared with KM .

  The effect on the product  P of  having a large or a small  substrate S is examined by considering the steady state value of I.

From (8.3) one has that the steady state value of I is

 

Isteady state = S ET/{S + (k2+k3)/k1)}                             

 

              =  S ET/{S + KM}                             .               (8.7)

 

Inserting the result in (8.5) gives the steady state rate of product growth

 

dP/dt  = k3 * S ET/{S + KM}                                          (8.8).

 

If S << KM  ,   

 

                      dP/dt  ≈ k3 * S ET/KM                      ,        (8.9)    

and the reaction is called substrate limited .                            

 

If S >> KM  ,

 

       dP/dt  ≈ k3  ET                                                            (8.10)

 

and the reaction is enzyme limited.                                      

 

 

 

 

Data for figure

         data et /.05/

      data k1,k2,k3/.5, 1.,10./

      data  int0, p0/0.,0./

Fig 17 .  S=0.1 KM  .

 

 

 

 

Data used  for fig 18.

     data et /.05/

      data k1,k2,k3/.5, 1.,10./

      data  int0, p0/0.,0./

Fig 18 .  S=5 KM .

 

The steady state product concentration is slightly  larger than  in Fig 18 and  takes a longer time to build up.

 

 

 

FORTRAN CODE

c Simon page 156 substrate and enzyme dependent reactions

c the vaules of k1,k2,k3 S , Et are arbitrary...They do not represent

c tyical value from physiology problems

      real k1,k2,k3,int0,int1

      data et /.05/

      data k1,k2,k3/.5, 1.,10./

      data  int0, p0/0.,0./

c S is substrate concentration

c low S is substrate limited reaction

      s=5.*(k2+k3)/k1

      Print*,'S=', s

      print*,'  '

      t2=1./k2

      t3=1./k3

      t1=1./(k1*s)

      tsmall= amin1(t1,t2,t3)

      tlarge=amax1(t1,t2,t3)

      dt=tsmall/500.

      tfinal=.09

      nstep=int(tfinal/dt)

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

      kount=kp

      print 100,0.,int0,p0

      do 10 i=1,nstep

      t=dt*float(i)

      a1=k1*s*(et-int0)-(k2+k3)*int0

      int1=int0+dt*a1

      p1=p0+dt*k3*int0

      if(i.eq.kount)then

      print 100,t,int1,p1

      kount=kount+kp

      endif

      int0=int1

      p1=p0

10    continue

100   format(1x,'t,I,P=',3(3x,e10.3))

      stop

      end

 

     

      END OF LECTURE