#file circadian_rhythm.m
#Goldbeter model of circadian rhythm in Drosophila
#reviewed in Goldbeter (1996) Biochemical oscillations and cellular
#rhythms, Cambridge University Press
#Figure 7.19

par vs=0.76
par vm=0.65
par vd=0.95
par ks=0.38
par kt1=1.9
par kt2=1.3
par v1=3.2
par v2=1.58
par v3=5
par v4=2.5
par k1=1
par k2=1
par k3=2
par k4=2
par ki=1
par km1=0.5
par kd=0.2
par n=4

#model equations
M' = vs/(1+(PN/ki)^n) - (vm*M)/(km1+M)
P0' = ks*M - (v1*P0)/(k1+P0) + (v2*P1)/(k2 + P1)
P1' = (v1*P0)/(k1+P0) - (v2*P1)/(k2 + P1) - (v3*P1)/(k3+P1) + (v4*P2)/(k4+P2)
P2' = (v3*P1)/(k3+P1) - (v4*P2)/(k4+P2)  - kt1*P2 + kt2*PN - (vd*P2)/(kd+P2)
PN' = kt1*P2 - kt2*PN

# initial condition
init M=1, P0=1, P1=0, P2=0, PN=0

# final time
@ total=90
# set the plotting window size:
@ xlo=0, xhi=90, ylo=0, yhi=4

# run the simulation in the GUI by selecting (I)nitialconds|(G)o
# plot other species by selecting (G)raphic stuff|(A)dd curve and 
# assigning the y-axis

#multiple trajectories can be generated by (G)raphic stuff|(F)reeze and
#running simulations from different initial conditions

#generate a phase plane by choosing (V)iewaxis|(2)D and assigning the axis
#the nullclines can be generated with (N)ullclies|(N)ew
#a direction field can be generated with (D)ir.field/flow|(D)irect field

done