# Stochastic version of Goodwin Oscillator
# V = Volume, M # molecules of mRNA, Po # molecules of protein, Pt # molecules of activated protein

init tr=0, M=10, MM=10, G=0

!azeroV=azero*V
!aoneV=aone*V
!binV= bin/V
!bMV= bM/V

p1 = (1-G)*azeroV + G*aoneV 
p2 = p1 + delta*M 
p3 = p2 + bMV*M*(M-1)
p4 = p3 + unbM*MM  
p5 = p4 + unbin*G
p6 = p5 + binV*(1-G)*MM

s2 = ran(1)*p6
z1 = (s2<p1)
z2 = (s2>=p1)&(s2<p2)
z3 = (s2>=p2)&(s2<p3)
z4 = (s2>=p3)&(s2<p4)
z5 = (s2>=p4)&(s2<p5)
z6 = (s2>=p5)&(s2<p6)

tr' = tr - log(ran(1))/p6
M' = max(0,M-z2+z1-2*z3+2*z4)
MM'= max(0,MM+z3-z4-z6+z5)
G' = max(0,G-z5+z6)


par azero=5.0, aone=220.0, delta=1.0, bM=1.0, unbM=100.0, bin=0.05, unbin=5.0
par V=1


aux MV=M/V
aux MMV=MM/V

@ meth=discrete, bound=100000000, total=500000, njmp=1
@ xp=tr, yp=M, xlo=0, xhi=1100, ylo=0, yhi=200
@ maxstor=40000000

done