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

init tr=0, M=5, Po=5, Pt=5, G=0

!transV=trans*V
!binV= bin/V


p1 = max(0, 1-G)*transV 
p2 = p1 + degm*M 
p3 = p2 + degpo*Po
p4 = p3 + degpt*Pt 
p5 = p4 + tl*M
p6 = p5 + conv*Po 
p7 = p6 + unbin*G
p8 = p7 + binV*(numsites-G)*Pt



s2 = ran(1)*p8
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)
z7 = (s2>=p6)&(s2<p7)
z8 = s2>=p7

tr' = tr - log(ran(1))/p8
M' = max(0,M-z2+z1)
Po' = max(0,Po-z3+z5-z6)
Pt' = max(0,Pt-z4+z6)
G' = max(0,G-z7+z8)


par trans=1.0, tl=1.0, degm=1.0, degpo=0.5, degpt=1.0, conv=0.5, bin=1000.0, unbin=10.0, numsites=4
par V=200


aux MV=M/V
aux PoV=Po/V
aux PtV=Pt/V

@ meth=discrete, bound=100000, total=40000, njmp=1
@ xp=tr, yp=M, xlo=0, xhi=100, ylo=0, yhi=75
@ maxstor=40000000

done