# bruss_gil.ode
# gillespie algorithm for brusselator
#
# here are the four reactions
# A -> U (c1)
# B+U -> V (c2)
# 2 U + V -> 3 U (c3)
# U -> * (c4)
#
par c1a=5000,c2b=50,c3=.00005,c4=5
init u=1000,v=2000
# compute the cumulative reaction rates
#
p1=c1a
p2=p1+c2b*u
p3=p2+c3*u*v*(u-1)/2
p4=p3+c4*u
# choose random # to determine which reaction
# note that p4 is the total of all the rates
s2=ran(1)*p4
# z1,z2,z3,z4 are either 0 or 1 according to whether the reaction
# took place or not
# if s2=p1)
z3=(s2=p2)
# finally, if s2 >= p3, then reaction 4
z4=(s2>=p3)
# time for next reaction
tr'=tr-log(ran(1))/p4
u'=max(0,u+z1-z2+z3-z4)
v'=max(0,v+z2-z3)
# this plots a phaseplane of the variables - it does 1000000 reactions
# and ploats out every thousandth
@ bound=100000000,meth=discrete,total=1000000,nout=1000
@ xp=u,yp=v
@ xlo=0,ylo=0,xhi=10000,yhi=10000
done