# file synchronized_repressilators.ode
# model of synchronized repressilator networks
# from Garcia-Ojalvo et al. (2004) PNAS 101 pp. 10955-10960
# problem 7.8.19

par alpha0=0 
par alpha=216  
par n=2
par beta=1

par kappa=20
par ks0=1
par eta=2
par ks1=0.01
par Q=0.9

#model equations
se=Q*((1/2)*(s1+s2))

ma1' = alpha0 + alpha/(1+pc1^n) - ma1
pa1' = beta*(ma1-pa1)
mb1' = alpha0 + alpha/(1+pa1^n) - mb1
pb1' = beta*(mb1-pb1)
mc1' = alpha0 + alpha/(1+pb1^n) + kappa*s1/(1+s1) - mc1
pc1' = beta*(mc1-pc1)
s1' = -ks0*s1+ks1*pa1-eta*(s1-se)

ma2' = alpha0 + alpha/(1+pc2^n) - ma2
pa2' = beta*(ma2-pa2)
mb2' = alpha0 + alpha/(1+pa2^n) - mb2
pb2' = beta*(mb2-pb2)
mc2' = alpha0 + alpha/(1+pb2^n) + kappa*s2/(1+s2) - mc2
pc2' = beta*(mc2-pc2)
s2' = -ks0*s2+ks1*pa2-eta*(s2-se)

# initial condition
init ma1=0, pa1=10, mb1=0, pb1=0, mc1=0, pc1=0, s1=0 
init ma2=0, pa2=0, mb2=0, pb2=10, mc1=0, pc1=0, s2=0 


# final time
@ total=100
# set the plotting window size:
@ xlo=0, xhi=100, ylo=0, yhi=100
@ bound=10000
# run the simulation in the GUI by selecting (I)nitialconds|(G)o
# plot other species by selecting (G)raphic stuff|(A)dd curve and 
# assigning the y-axis

done