#file stochastic_gene_expression.ode
#stochastic implementation (Gillespie SSA) of constitutive gene expression
#Figure 7.46

#set propensity parameters

#burst size
par b=1

par kr=10, kp=6, gr=1, gp=1

#set initial condition for molecular counts M, P , and for elapsed time tr.
#
init M=0, P=0, tr=0
#
#set reaction propensity
#
a1=kr
a2=kp*M
a3=gr*M
a4=gp*P

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)

M' = M+b*z0 - z2
P' = P+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=1000
#
# 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=140
#
done
#