#file stochastic_brusselator.ode
#stochastic implementation (Gillespie SSA) of the brusselator model
#Figure 7.47

#set propensity parameters

par k1=5000, k2=50, k3=0.00005, k4=5

#set initial condition for molecular counts X, Y , and for elapsed time tr.
#
init X=1000, Y=2000, tr=0
#
#set reaction propensity
#

a1=k1
a2=k2*X
a3=k3*Y*X*(X-1)/2
a4=k4*X

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-z1+z2-z3
Y' = Y+z1-z2

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


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