# file stochastic_toggle_switch.ode
# stochastic simulation of genetic toggle switch
# problem 7.8.26

#set propensity parameters
par d=1
par a=5
par Km=1
par beta=4

#set initial condition for molecular counts M, P , and for elapsed time tr.
#
init X=0, Y=5, tr=0
#
#set reaction propensity
#
a1=a/(Km+Y)
a2=a/(Km+X)
a3=d*X
a4=d*Y

asum=a1+a2+a3+a4

p1=a1/asum
p2=(a1+a2)/asum
p3=(a1+a2+a3)/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<=1)

X' = X+z0 - z2
Y' = Y+z1-z3

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

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