# file stochastic_vilar_oscillator.ode
# stochastic simulation of genetic oscillator
# from Vilar et al. (2002) PNAS 99 5988-5992.
# problem 7.8.27

#set propensity parameters

ga=250
ba=5
Ka=0.5 
alpha0=0.1
da=1 
kc=200 
gr=50 
br=10
Kr=1
dr=0.1

#set initial condition for molecular counts A, R, C and for elapsed time tr.
#
init A=0, R=0, C=0, tr=0
#
#set reaction propensity
#
a1=(ga/ba)*(alpha0+(A/Ka))/(1+(A/Ka))
a2=da*A
a3=kc*A*R
a4=(gr/br)*((A/Kr))/(1+(A/Kr))
a5=dr*R
a6=da*C

asum=a1+a2+a3+a4+a5+a6

p1=a1/asum
p2=(a1+a2)/asum
p3=(a1+a2+a3)/asum
p4=(a1+a2+a3+a4)/asum
p5=(a1+a2+a3+a4+a5)/asum


#set the update rules for molecule abundance
mu=ran(1)

z0=(0<=mu)&(mu<p1)
z1=(p1<=mu)&(mu<p2)
z2=(p2<=mu)&(mu<p3)
z3=(p3<=mu)&(mu<p4)
z4=(p4<=mu)&(mu<p5)
z5=(p5<=mu)&(mu<=1)

A' = A + ba*z0 - z1 -z2
R' = R + br*z3 -z4 -z2+z5
C' = C + z2 -z5

# set the update rule for tr
#
tr'=tr+log(1/ran(1))/asum
#

# set the total number of steps to be taken
@ total=20000
@ maxstor=1000000
@ bound=1000
#
# set the model as discete-time (update rule-based)
@ meth=discrete
#
# set the plotting default to be A against tr
@ xp=tr,yp=A, xlo=0,xhi=70,ylo=0,yhi=200
#
done
#