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:
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 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 DEs 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