#file genetic_relaxation_oscillator.ode
#Hasty synthetic oscillator model
#from Hasty et al. (2002) Phys. Rev. Lett. 88 148101
#Figure 7.23
#NOTE: typos in the published book: 
#page 263, caption of figure 7.23, the value for gamma_y should be 0.012, not 0.06.
#page 303 Problem 7.8.14, the value of gamma_y should be  0.15, (not 1 as written).

par alpha=11
par sigma=2
par gammax=0.2
par gammay=0.012
par ay=0.2

#model equations
x' = (1+x^2+alpha*sigma*x^4)/((1+x^2+sigma*x^4)*(1+y^4)) - gammax*x
y' = ay*((1+x^2+alpha*sigma*x^4)/((1+x^2+sigma*x^4)*(1+y^4))) - gammay*y


# initial condition
init x=0.3963, y=2.3346

# final time
@ total=300
# set the plotting window size:
@ xlo=0, xhi=300, ylo=0, yhi=3
@ maxstor=100000

# run the simulation in the GUI by selecting (I)nitialconds|(G)o
# plot species y by selecting (G)raphic stuff|(A)dd curve and 
# assigning the y-axis

#generate a phase plane by choosing (V)iewaxis|(2)D and assigning the axis
#the nullclines can be generated with (N)ullclies|(N)ew
#a direction field can be generated with (D)ir.field/flow|(D)irect field

done